libMesh
nemesis_io_helper.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 // C++ headers
20 #include <iomanip>
21 #include <set>
22 #include <sstream>
23 
24 // Libmesh headers
25 #include "libmesh/nemesis_io_helper.h"
26 
27 #include "libmesh/boundary_info.h"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/elem.h"
30 #include "libmesh/equation_systems.h"
31 #include "libmesh/fe_interface.h"
32 #include "libmesh/int_range.h"
33 #include "libmesh/mesh_base.h"
34 #include "libmesh/node.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/parallel.h"
37 #include "libmesh/utility.h"
38 
39 #if defined(LIBMESH_HAVE_NEMESIS_API) && defined(LIBMESH_HAVE_EXODUS_API)
40 
41 #include <libmesh/ignore_warnings.h>
42 namespace exII {
43 extern "C" {
44 #include "exodusII.h" // defines MAX_LINE_LENGTH, MAX_STR_LENGTH used later
45 }
46 }
47 #include <libmesh/restore_warnings.h>
48 
49 namespace libMesh
50 {
51 
52 
53 // Initialize the various integer members to zero. We can check
54 // these later to see if they've been properly initialized...
55 // The parent ExodusII_IO_Helper is created with the run_only_on_proc0
56 // flag set to false, so that we can make use of its functionality
57 // on multiple processors.
59  bool verbose_in, bool single_precision) :
60  ExodusII_IO_Helper(parent, verbose_in, /*run_only_on_proc0=*/false, /*single_precision=*/single_precision),
61  nemesis_err_flag(0),
62  num_nodes_global(0),
63  num_elems_global(0),
64  num_elem_blks_global(0),
65  num_node_sets_global(0),
66  num_side_sets_global(0),
67  num_proc(0),
68  num_proc_in_file(0),
69  ftype('\0'),
70  num_internal_nodes(0),
71  num_border_nodes(0),
72  num_external_nodes(0),
73  num_internal_elems(0),
74  num_border_elems(0),
75  num_node_cmaps(0),
76  num_elem_cmaps(0),
77  write_complex_abs(true)
78 {
79  // Warn about using untested code!
80  libmesh_experimental();
81 }
82 
83 
85 {
86  // Our destructor is called from Nemesis_IO. We close the Exodus file here since we have
87  // responsibility for managing the file's lifetime. Only call ex_update() if the file was
88  // opened for writing!
89  if (this->opened_for_writing)
90  {
91  this->ex_err = exII::ex_update(this->ex_id);
92  EX_EXCEPTIONLESS_CHECK_ERR(ex_err, "Error flushing buffers to file.");
93  }
94  this->close();
95 }
96 
97 
98 
100 {
101  libmesh_assert_less (id, nodeset_ids.size());
102  libmesh_assert_less (id, num_nodes_per_set.size());
103  libmesh_assert_less (id, num_node_df_per_set.size());
104 
105  ex_err = exII::ex_get_set_param(ex_id,
106  exII::EX_NODE_SET,
107  nodeset_ids[id],
108  &num_nodes_per_set[id],
109  &num_node_df_per_set[id]);
110  EX_CHECK_ERR(ex_err, "Error retrieving nodeset parameters.");
111  message("Parameters retrieved successfully for nodeset: ", id);
112 
113  node_list.resize(num_nodes_per_set[id]);
114 
115  // Don't call ex_get_set unless there are actually nodes there to get.
116  // Exodus prints an annoying warning message in DEBUG mode otherwise...
117  if (num_nodes_per_set[id] > 0)
118  {
119  ex_err = exII::ex_get_set(ex_id,
120  exII::EX_NODE_SET,
121  nodeset_ids[id],
122  node_list.data(),
123  nullptr); // set_extra_list, ignored for node sets
124 
125  EX_CHECK_ERR(ex_err, "Error retrieving nodeset data.");
126  message("Data retrieved successfully for nodeset: ", id);
127  }
128 }
129 
130 
131 
133 {
135  Nemesis::ne_get_init_global(ex_id,
141  EX_CHECK_ERR(nemesis_err_flag, "Error reading initial global data!");
142 
143  if (verbose)
144  {
145  libMesh::out << "[" << this->processor_id() << "] " << "num_nodes_global=" << num_nodes_global << std::endl;
146  libMesh::out << "[" << this->processor_id() << "] " << "num_elems_global=" << num_elems_global << std::endl;
147  libMesh::out << "[" << this->processor_id() << "] " << "num_elem_blks_global=" << num_elem_blks_global << std::endl;
148  libMesh::out << "[" << this->processor_id() << "] " << "num_node_sets_global=" << num_node_sets_global << std::endl;
149  libMesh::out << "[" << this->processor_id() << "] " << "num_side_sets_global=" << num_side_sets_global << std::endl;
150  }
151 }
152 
153 
154 
156 {
157  if (num_side_sets_global > 0)
158  {
161 
162  // df = "distribution factor", not really sure what that is. I don't yet have a file
163  // which has distribution factors so I guess we'll worry about it later...
165 
167  Nemesis::ne_get_ss_param_global(ex_id,
168  global_sideset_ids.data(),
169  num_global_side_counts.data(),
171  EX_CHECK_ERR(nemesis_err_flag, "Error reading global sideset parameters!");
172 
173  if (verbose)
174  {
175  libMesh::out << "[" << this->processor_id() << "] " << "Global Sideset IDs, Side Counts, and DF counts:" << std::endl;
176  for (auto bn : index_range(global_sideset_ids))
177  {
178  libMesh::out << " [" << this->processor_id() << "] "
179  << "global_sideset_ids["<<bn<<"]=" << global_sideset_ids[bn]
180  << ", num_global_side_counts["<<bn<<"]=" << num_global_side_counts[bn]
181  << ", num_global_side_df_counts["<<bn<<"]=" << num_global_side_df_counts[bn]
182  << std::endl;
183  }
184  }
185  }
186 }
187 
188 
189 
190 
192 {
193  if (num_node_sets_global > 0)
194  {
198 
200  Nemesis::ne_get_ns_param_global(ex_id,
201  global_nodeset_ids.data(),
202  num_global_node_counts.data(),
204  EX_CHECK_ERR(nemesis_err_flag, "Error reading global nodeset parameters!");
205 
206  if (verbose)
207  {
208  libMesh::out << "[" << this->processor_id() << "] " << "Global Nodeset IDs, Node Counts, and DF counts:" << std::endl;
209  for (auto bn : index_range(global_nodeset_ids))
210  {
211  libMesh::out << " [" << this->processor_id() << "] "
212  << "global_nodeset_ids["<<bn<<"]=" << global_nodeset_ids[bn]
213  << ", num_global_node_counts["<<bn<<"]=" << num_global_node_counts[bn]
214  << ", num_global_node_df_counts["<<bn<<"]=" << num_global_node_df_counts[bn]
215  << std::endl;
216  }
217  }
218  }
219 }
220 
221 
222 
224 {
227 
228  if (num_elem_blks_global > 0)
229  {
231  Nemesis::ne_get_eb_info_global(ex_id,
232  global_elem_blk_ids.data(),
233  global_elem_blk_cnts.data());
234  EX_CHECK_ERR(nemesis_err_flag, "Error reading global element block info!");
235  }
236 
237  if (verbose)
238  {
239  libMesh::out << "[" << this->processor_id() << "] " << "Global Element Block IDs and Counts:" << std::endl;
240  for (auto bn : index_range(global_elem_blk_ids))
241  {
242  libMesh::out << " [" << this->processor_id() << "] "
243  << "global_elem_blk_ids[" << bn << "]=" << global_elem_blk_ids[bn]
244  << ", global_elem_blk_cnts[" << bn << "]=" << global_elem_blk_cnts[bn]
245  << std::endl;
246  }
247  }
248 }
249 
250 
251 
253 {
255  Nemesis::ne_get_init_info(ex_id,
256  &num_proc,
258  &ftype);
259  EX_CHECK_ERR(nemesis_err_flag, "Error reading initial info!");
260 
261  if (verbose)
262  {
263  libMesh::out << "[" << this->processor_id() << "] " << "num_proc=" << num_proc << std::endl;
264  libMesh::out << "[" << this->processor_id() << "] " << "num_proc_in_file=" << num_proc_in_file << std::endl;
265  libMesh::out << "[" << this->processor_id() << "] " << "ftype=" << ftype << std::endl;
266  }
267 }
268 
269 
270 
272 {
274  Nemesis::ne_get_loadbal_param(ex_id,
282  this->processor_id() // The ID of the processor for which info is to be read
283  );
284  EX_CHECK_ERR(nemesis_err_flag, "Error reading load balance parameters!");
285 
286 
287  if (verbose)
288  {
289  libMesh::out << "[" << this->processor_id() << "] " << "num_internal_nodes=" << num_internal_nodes << std::endl;
290  libMesh::out << "[" << this->processor_id() << "] " << "num_border_nodes=" << num_border_nodes << std::endl;
291  libMesh::out << "[" << this->processor_id() << "] " << "num_external_nodes=" << num_external_nodes << std::endl;
292  libMesh::out << "[" << this->processor_id() << "] " << "num_internal_elems=" << num_internal_elems << std::endl;
293  libMesh::out << "[" << this->processor_id() << "] " << "num_border_elems=" << num_border_elems << std::endl;
294  libMesh::out << "[" << this->processor_id() << "] " << "num_node_cmaps=" << num_node_cmaps << std::endl;
295  libMesh::out << "[" << this->processor_id() << "] " << "num_elem_cmaps=" << num_elem_cmaps << std::endl;
296  }
297 }
298 
299 
300 
302 {
304  elem_mapb.resize(num_border_elems);
305 
307  Nemesis::ne_get_elem_map(ex_id,
308  elem_mapi.empty() ? nullptr : elem_mapi.data(),
309  elem_mapb.empty() ? nullptr : elem_mapb.data(),
310  this->processor_id()
311  );
312  EX_CHECK_ERR(nemesis_err_flag, "Error reading element maps!");
313 
314 
315  if (verbose)
316  {
317  libMesh::out << "[" << this->processor_id() << "] elem_mapi[i] = ";
318  for (auto i : make_range(num_internal_elems-1))
319  libMesh::out << elem_mapi[i] << ", ";
320  libMesh::out << "... " << elem_mapi.back() << std::endl;
321 
322  libMesh::out << "[" << this->processor_id() << "] elem_mapb[i] = ";
323  for (auto i : make_range(std::min(10, num_border_elems-1)))
324  libMesh::out << elem_mapb[i] << ", ";
325  libMesh::out << "... " << elem_mapb.back() << std::endl;
326  }
327 }
328 
329 
330 
331 
333 {
335  node_mapb.resize(num_border_nodes);
337 
339  Nemesis::ne_get_node_map(ex_id,
340  node_mapi.empty() ? nullptr : node_mapi.data(),
341  node_mapb.empty() ? nullptr : node_mapb.data(),
342  node_mape.empty() ? nullptr : node_mape.data(),
343  this->processor_id()
344  );
345  EX_CHECK_ERR(nemesis_err_flag, "Error reading node maps!");
346 
347  if (verbose)
348  {
349  // Remark: The Exodus/Nemesis node numbering is always (?) 1-based! This means the first interior node id will
350  // always be == 1.
351  libMesh::out << "[" << this->processor_id() << "] " << "first interior node id=" << node_mapi[0] << std::endl;
352  libMesh::out << "[" << this->processor_id() << "] " << "last interior node id=" << node_mapi.back() << std::endl;
353 
354  libMesh::out << "[" << this->processor_id() << "] " << "first boundary node id=" << node_mapb[0] << std::endl;
355  libMesh::out << "[" << this->processor_id() << "] " << "last boundary node id=" << node_mapb.back() << std::endl;
356 
357  // The number of external nodes is sometimes zero, don't try to access
358  // node_mape.back() in this case!
359  if (num_external_nodes > 0)
360  {
361  libMesh::out << "[" << this->processor_id() << "] " << "first external node id=" << node_mape[0] << std::endl;
362  libMesh::out << "[" << this->processor_id() << "] " << "last external node id=" << node_mape.back() << std::endl;
363  }
364  }
365 }
366 
367 
368 
370 {
375 
377  Nemesis::ne_get_cmap_params(ex_id,
378  node_cmap_ids.empty() ? nullptr : node_cmap_ids.data(),
379  node_cmap_node_cnts.empty() ? nullptr : node_cmap_node_cnts.data(),
380  elem_cmap_ids.empty() ? nullptr : elem_cmap_ids.data(),
381  elem_cmap_elem_cnts.empty() ? nullptr : elem_cmap_elem_cnts.data(),
382  this->processor_id());
383  EX_CHECK_ERR(nemesis_err_flag, "Error reading cmap parameters!");
384 
385 
386  if (verbose)
387  {
388  libMesh::out << "[" << this->processor_id() << "] ";
389  for (auto i : index_range(node_cmap_ids))
390  libMesh::out << "node_cmap_ids[" << i << "]=" << node_cmap_ids[i] << " ";
391  libMesh::out << std::endl;
392 
393  libMesh::out << "[" << this->processor_id() << "] ";
394  for (auto i : index_range(node_cmap_node_cnts))
395  libMesh::out << "node_cmap_node_cnts[" << i << "]=" << node_cmap_node_cnts[i] << " ";
396  libMesh::out << std::endl;
397 
398  libMesh::out << "[" << this->processor_id() << "] ";
399  for (auto i : index_range(elem_cmap_ids))
400  libMesh::out << "elem_cmap_ids[" << i << "]=" << elem_cmap_ids[i] << " ";
401  libMesh::out << std::endl;
402 
403  libMesh::out << "[" << this->processor_id() << "] ";
404  for (auto i : index_range(elem_cmap_elem_cnts))
405  libMesh::out << "elem_cmap_elem_cnts[" << i << "]=" << elem_cmap_elem_cnts[i] << " ";
406  libMesh::out << std::endl;
407  }
408 }
409 
410 
411 
413 {
416 
417  for (auto i : index_range(node_cmap_node_ids))
418  {
421 
422  // Don't call ne_get_node_cmap() if there is nothing there to
423  // get, Nemesis throws an error in this case.
424  if (node_cmap_node_cnts[i] > 0)
425  {
427  Nemesis::ne_get_node_cmap(ex_id,
428  node_cmap_ids[i],
429  node_cmap_node_ids[i].data(),
430  node_cmap_proc_ids[i].data(),
431  this->processor_id());
432  EX_CHECK_ERR(nemesis_err_flag, "Error reading node cmap node and processor ids!");
433  }
434 
435  if (verbose)
436  {
437  libMesh::out << "[" << this->processor_id() << "] node_cmap_node_ids[" << i << "]=";
438  for (const auto & dof : node_cmap_node_ids[i])
439  libMesh::out << dof << " ";
440  libMesh::out << std::endl;
441 
442  // This is basically a vector, all entries of which are = node_cmap_ids[i]
443  // Not sure if it's always guaranteed to be that or what...
444  libMesh::out << "[" << this->processor_id() << "] node_cmap_proc_ids[" << i << "]=";
445  for (const auto & dof : node_cmap_proc_ids[i])
446  libMesh::out << dof << " ";
447  libMesh::out << std::endl;
448  }
449  }
450 }
451 
452 
453 
455 {
459 
460  for (auto i : index_range(elem_cmap_elem_ids))
461  {
465 
466  if (elem_cmap_elem_cnts[i] > 0)
467  {
469  Nemesis::ne_get_elem_cmap(ex_id,
470  elem_cmap_ids[i],
471  elem_cmap_elem_ids[i].data(),
472  elem_cmap_side_ids[i].data(),
473  elem_cmap_proc_ids[i].data(),
474  this->processor_id());
475  EX_CHECK_ERR(nemesis_err_flag, "Error reading elem cmap elem, side, and processor ids!");
476  }
477 
478  if (verbose)
479  {
480  libMesh::out << "[" << this->processor_id() << "] elem_cmap_elem_ids[" << i << "]=";
481  for (const auto & dof : elem_cmap_elem_ids[i])
482  libMesh::out << dof << " ";
483  libMesh::out << std::endl;
484 
485  // These must be the (local) side IDs (in the ExodusII face numbering scheme)
486  // of the sides shared across processors.
487  libMesh::out << "[" << this->processor_id() << "] elem_cmap_side_ids[" << i << "]=";
488  for (const auto & dof : elem_cmap_side_ids[i])
489  libMesh::out << dof << " ";
490  libMesh::out << std::endl;
491 
492  // This is basically a vector, all entries of which are = elem_cmap_ids[i]
493  // Not sure if it's always guaranteed to be that or what...
494  libMesh::out << "[" << this->processor_id() << "] elem_cmap_proc_ids[" << i << "]=";
495  for (const auto & dof : elem_cmap_proc_ids[i])
496  libMesh::out << dof << " ";
497  libMesh::out << std::endl;
498  }
499  }
500 }
501 
502 
503 
504 
505 void Nemesis_IO_Helper::put_init_info(unsigned num_proc_in,
506  unsigned num_proc_in_file_in,
507  const char * ftype_in)
508 {
510  Nemesis::ne_put_init_info(ex_id,
511  num_proc_in,
512  num_proc_in_file_in,
513  const_cast<char *>(ftype_in));
514 
515  EX_CHECK_ERR(nemesis_err_flag, "Error writing initial information!");
516 }
517 
518 
519 
520 
522  dof_id_type num_elems_global_in,
523  unsigned num_elem_blks_global_in,
524  unsigned num_node_sets_global_in,
525  unsigned num_side_sets_global_in)
526 {
528  Nemesis::ne_put_init_global(ex_id,
529  num_nodes_global_in,
530  num_elems_global_in,
531  num_elem_blks_global_in,
532  num_node_sets_global_in,
533  num_side_sets_global_in);
534 
535  EX_CHECK_ERR(nemesis_err_flag, "Error writing initial global data!");
536 }
537 
538 
539 
540 void Nemesis_IO_Helper::put_eb_info_global(std::vector<int> & global_elem_blk_ids_in,
541  std::vector<int> & global_elem_blk_cnts_in)
542 {
544  Nemesis::ne_put_eb_info_global(ex_id,
545  global_elem_blk_ids_in.data(),
546  global_elem_blk_cnts_in.data());
547 
548  EX_CHECK_ERR(nemesis_err_flag, "Error writing global element block information!");
549 }
550 
551 
552 
553 
554 void Nemesis_IO_Helper::put_ns_param_global(std::vector<int> & global_nodeset_ids_in,
555  std::vector<int> & num_global_node_counts_in,
556  std::vector<int> & num_global_node_df_counts_in)
557 {
558  // Only add nodesets if there are some
559  if (global_nodeset_ids.size())
560  {
562  Nemesis::ne_put_ns_param_global(ex_id,
563  global_nodeset_ids_in.data(),
564  num_global_node_counts_in.data(),
565  num_global_node_df_counts_in.data());
566  }
567 
568  EX_CHECK_ERR(nemesis_err_flag, "Error writing global nodeset parameters!");
569 }
570 
571 
572 
573 
574 void Nemesis_IO_Helper::put_ss_param_global(std::vector<int> & global_sideset_ids_in,
575  std::vector<int> & num_global_side_counts_in,
576  std::vector<int> & num_global_side_df_counts_in)
577 {
578  // Only add sidesets if there are some
579  if (global_sideset_ids.size())
580  {
582  Nemesis::ne_put_ss_param_global(ex_id,
583  global_sideset_ids_in.data(),
584  num_global_side_counts_in.data(),
585  num_global_side_df_counts_in.data());
586  }
587 
588  EX_CHECK_ERR(nemesis_err_flag, "Error writing global sideset parameters!");
589 }
590 
591 
592 
593 
594 void Nemesis_IO_Helper::put_loadbal_param(unsigned num_internal_nodes_in,
595  unsigned num_border_nodes_in,
596  unsigned num_external_nodes_in,
597  unsigned num_internal_elems_in,
598  unsigned num_border_elems_in,
599  unsigned num_node_cmaps_in,
600  unsigned num_elem_cmaps_in)
601 {
603  Nemesis::ne_put_loadbal_param(ex_id,
604  num_internal_nodes_in,
605  num_border_nodes_in,
606  num_external_nodes_in,
607  num_internal_elems_in,
608  num_border_elems_in,
609  num_node_cmaps_in,
610  num_elem_cmaps_in,
611  this->processor_id());
612 
613  EX_CHECK_ERR(nemesis_err_flag, "Error writing loadbal parameters!");
614 }
615 
616 
617 
618 
619 
620 void Nemesis_IO_Helper::put_cmap_params(std::vector<int> & node_cmap_ids_in,
621  std::vector<int> & node_cmap_node_cnts_in,
622  std::vector<int> & elem_cmap_ids_in,
623  std::vector<int> & elem_cmap_elem_cnts_in)
624 {
626  Nemesis::ne_put_cmap_params(ex_id,
627  node_cmap_ids_in.empty() ? nullptr : node_cmap_ids_in.data(),
628  node_cmap_node_cnts_in.empty() ? nullptr : node_cmap_node_cnts_in.data(),
629  elem_cmap_ids_in.empty() ? nullptr : elem_cmap_ids_in.data(),
630  elem_cmap_elem_cnts_in.empty() ? nullptr : elem_cmap_elem_cnts_in.data(),
631  this->processor_id());
632 
633  EX_CHECK_ERR(nemesis_err_flag, "Error writing cmap parameters!");
634 }
635 
636 
637 
638 
639 void Nemesis_IO_Helper::put_node_cmap(std::vector<std::vector<int>> & node_cmap_node_ids_in,
640  std::vector<std::vector<int>> & node_cmap_proc_ids_in)
641 {
642  // Print to screen what we are about to print to Nemesis file
643  if (verbose)
644  {
645  for (auto i : index_range(node_cmap_node_ids_in))
646  {
647  libMesh::out << "[" << this->processor_id() << "] put_node_cmap() : nodes communicated to proc "
648  << this->node_cmap_ids[i]
649  << " = ";
650  for (const auto & node_id : node_cmap_node_ids_in[i])
651  libMesh::out << node_id << " ";
652  libMesh::out << std::endl;
653  }
654 
655  for (auto i : index_range(node_cmap_node_ids_in))
656  {
657  libMesh::out << "[" << this->processor_id() << "] put_node_cmap() : processor IDs = ";
658  for (const auto & proc_id : node_cmap_proc_ids_in[i])
659  libMesh::out << proc_id << " ";
660  libMesh::out << std::endl;
661  }
662  }
663 
664  for (auto i : index_range(node_cmap_node_ids_in))
665  {
666  int * node_ids_ptr = node_cmap_node_ids_in[i].empty() ?
667  nullptr : node_cmap_node_ids_in[i].data();
668  int * proc_ids_ptr = node_cmap_proc_ids_in[i].empty() ?
669  nullptr : node_cmap_proc_ids_in[i].data();
670 
672  Nemesis::ne_put_node_cmap(ex_id, this->node_cmap_ids[i],
673  node_ids_ptr, proc_ids_ptr,
674  this->processor_id());
675 
676  EX_CHECK_ERR(nemesis_err_flag, "Error writing node communication map to file!");
677  }
678 }
679 
680 
681 
682 
683 void Nemesis_IO_Helper::put_node_map(std::vector<int> & node_mapi_in,
684  std::vector<int> & node_mapb_in,
685  std::vector<int> & node_mape_in)
686 {
688  Nemesis::ne_put_node_map(ex_id,
689  node_mapi_in.empty() ? nullptr : node_mapi_in.data(),
690  node_mapb_in.empty() ? nullptr : node_mapb_in.data(),
691  node_mape_in.empty() ? nullptr : node_mape_in.data(),
692  this->processor_id());
693 
694  EX_CHECK_ERR(nemesis_err_flag, "Error writing Nemesis internal and border node maps to file!");
695 }
696 
697 
698 
699 
700 void Nemesis_IO_Helper::put_elem_cmap(std::vector<std::vector<int>> & elem_cmap_elem_ids_in,
701  std::vector<std::vector<int>> & elem_cmap_side_ids_in,
702  std::vector<std::vector<int>> & elem_cmap_proc_ids_in)
703 {
704  for (auto i : index_range(elem_cmap_ids))
705  {
707  Nemesis::ne_put_elem_cmap(ex_id,
708  this->elem_cmap_ids[i],
709  elem_cmap_elem_ids_in[i].data(),
710  elem_cmap_side_ids_in[i].data(),
711  elem_cmap_proc_ids_in[i].data(),
712  this->processor_id());
713 
714  EX_CHECK_ERR(nemesis_err_flag, "Error writing elem communication map to file!");
715  }
716 }
717 
718 
719 
720 
721 void Nemesis_IO_Helper::put_elem_map(std::vector<int> & elem_mapi_in,
722  std::vector<int> & elem_mapb_in)
723 {
725  Nemesis::ne_put_elem_map(ex_id,
726  elem_mapi_in.empty() ? nullptr : elem_mapi_in.data(),
727  elem_mapb_in.empty() ? nullptr : elem_mapb_in.data(),
728  this->processor_id());
729 
730  EX_CHECK_ERR(nemesis_err_flag, "Error writing Nemesis internal and border element maps to file!");
731 }
732 
733 
734 
735 void Nemesis_IO_Helper::initialize(std::string title_in, const MeshBase & mesh, bool /*use_discontinuous*/)
736 {
737  // Make sure that the reference passed in is really a DistributedMesh
738  // const DistributedMesh & pmesh = cast_ref<const DistributedMesh &>(mesh);
739  const MeshBase & pmesh = mesh;
740 
741  // If _write_as_dimension is nonzero, use it to set num_dim later in the Exodus file.
745  num_dim = mesh.mesh_dimension();
746  else
747  num_dim = mesh.spatial_dimension();
748 
749  // According to Nemesis documentation, first call when writing should be to
750  // ne_put_init_info(). Our reader doesn't actually call this, but we should
751  // strive to be as close to a normal nemesis file as possible...
752  this->put_init_info(this->n_processors(), 1, "p");
753 
754 
755  // Gather global "initial" information for Nemesis. This consists of
756  // three parts labeled I, II, and III below...
757 
758  //
759  // I.) Need to compute the number of global element blocks. To be consistent with
760  // Exodus, we also incorrectly associate the number of element blocks with the
761  // number of libmesh subdomains...
762  //
763  this->compute_num_global_elem_blocks(pmesh);
764 
765  //
766  // II.) Determine the global number of nodesets by communication.
767  // This code relies on BoundaryInfo storing side and node
768  // boundary IDs separately at the time they are added to the
769  // BoundaryInfo object.
770  //
771  this->compute_num_global_nodesets(pmesh);
772 
773  //
774  // III.) Need to compute the global number of sidesets by communication:
775  // This code relies on BoundaryInfo storing side and node
776  // boundary IDs separately at the time they are added to the
777  // BoundaryInfo object.
778  //
779  this->compute_num_global_sidesets(pmesh);
780 
781  // Now write the global data obtained in steps I, II, and III to the Nemesis file
782  this->put_init_global(pmesh.parallel_n_nodes(),
783  pmesh.parallel_n_elem(),
784  this->num_elem_blks_global, /* I. */
785  this->num_node_sets_global, /* II. */
786  this->num_side_sets_global /* III. */
787  );
788 
789  // Next, we'll write global element block information to the file. This was already
790  // gathered in step I. above
792  this->global_elem_blk_cnts);
793 
794 
795  // Next, write global nodeset information to the file. This was already gathered in
796  // step II. above.
797  this->num_global_node_df_counts.clear();
798  this->num_global_node_df_counts.resize(this->global_nodeset_ids.size()); // distribution factors all zero...
802 
803 
804  // Next, write global sideset information to the file. This was already gathered in
805  // step III. above.
806  this->num_global_side_df_counts.clear();
807  this->num_global_side_df_counts.resize(this->global_sideset_ids.size()); // distribution factors all zero...
811 
812 
813  // Before we go any further we need to derive consistent node and
814  // element numbering schemes for all local elems and nodes connected
815  // to local elements.
816  //
817  // Must be called *after* the local_subdomain_counts map has been constructed
818  // by the compute_num_global_elem_blocks() function!
819  this->build_element_and_node_maps(pmesh);
820 
821  // Next step is to write "load balance" parameters. Several things need to
822  // be computed first though...
823 
824  // First we'll collect IDs of border nodes.
825  this->compute_border_node_ids(pmesh);
826 
827  // Next we'll collect numbers of internal and border elements, and internal nodes.
828  // Note: "A border node does not a border element make...", that is, just because one
829  // of an element's nodes has been identified as a border node, the element is not
830  // necessarily a border element. It must have a side on the boundary between processors,
831  // i.e. have a face neighbor with a different processor id...
833 
834  // Finally we are ready to write the loadbal information to the file
836  this->num_border_nodes,
837  this->num_external_nodes,
838  this->num_internal_elems,
839  this->num_border_elems,
840  this->num_node_cmaps,
841  this->num_elem_cmaps);
842 
843 
844  // Now we need to compute the "communication map" parameters. These are basically
845  // lists of nodes and elements which need to be communicated between different processors
846  // when the mesh file is read back in.
848 
849  // Write communication map parameters to file.
850  this->put_cmap_params(this->node_cmap_ids,
851  this->node_cmap_node_cnts,
852  this->elem_cmap_ids,
853  this->elem_cmap_elem_cnts);
854 
855  // Ready the node communication maps. The node IDs which
856  // are communicated are the ones currently stored in
857  // proc_nodes_touched_intersections.
859 
860  // Write the packed node communication vectors to file.
861  this->put_node_cmap(this->node_cmap_node_ids,
862  this->node_cmap_proc_ids);
863 
864  // Ready the node maps. These have nothing to do with communication, they map
865  // the nodes to internal, border, and external nodes in the file.
866  this->compute_node_maps();
867 
868  // Call the Nemesis API to write the node maps to file.
869  this->put_node_map(this->node_mapi,
870  this->node_mapb,
871  this->node_mape);
872 
873  // Ready the element communication maps. This includes border
874  // element IDs, sides which are on the border, and the processors to which
875  // they are to be communicated...
877 
878  // Call the Nemesis API to write the packed element communication maps vectors to file
879  this->put_elem_cmap(this->elem_cmap_elem_ids,
880  this->elem_cmap_side_ids,
881  this->elem_cmap_proc_ids);
882 
883  // Ready the Nemesis element maps (internal and border) for writing to file.
884  this->compute_element_maps();
885 
886  // Call the Nemesis API to write the internal and border element IDs.
887  this->put_elem_map(this->elem_mapi,
888  this->elem_mapb);
889 
890  // Now write Exodus-specific initialization information, some of which is
891  // different when you are using Nemesis.
892  this->write_exodus_initialization_info(pmesh, title_in);
893 } // end initialize()
894 
895 
896 
897 
898 
899 
901  const std::string & title_in)
902 {
903  this->num_elem = static_cast<unsigned int>(std::distance (pmesh.active_local_elements_begin(),
904  pmesh.active_local_elements_end()));
905 
906  // Exodus will also use *global* number of side and node sets,
907  // though it will not write out entries for all of them...
908  this->num_side_sets =
909  cast_int<int>(this->global_sideset_ids.size());
910  this->num_node_sets =
911  cast_int<int>(this->global_nodeset_ids.size());
912 
913  // We need to write the global number of blocks, even though this processor might not have
914  // elements in some of them!
915  this->num_elem_blk = this->num_elem_blks_global;
916 
917  ex_err = exII::ex_put_init(ex_id,
918  title_in.c_str(),
919  this->num_dim,
920  this->num_nodes,
921  this->num_elem,
922  this->num_elem_blk,
923  this->num_node_sets,
924  this->num_side_sets);
925 
926  EX_CHECK_ERR(ex_err, "Error initializing new Nemesis file.");
927 }
928 
929 
930 
931 
932 
934 {
935  // Make sure we don't have any leftover info
936  this->elem_mapi.clear();
937  this->elem_mapb.clear();
938 
939  // Copy set contents into vectors
940  this->elem_mapi.resize(this->internal_elem_ids.size());
941  this->elem_mapb.resize(this->border_elem_ids.size());
942 
943  {
944  unsigned cnt = 0;
945  for (const auto & id : this->internal_elem_ids)
946  this->elem_mapi[cnt++] = libmesh_map_find(libmesh_elem_num_to_exodus, id);
947  }
948 
949  {
950  unsigned cnt = 0;
951  for (const auto & id : this->border_elem_ids)
952  this->elem_mapb[cnt++] = libmesh_map_find(libmesh_elem_num_to_exodus, id);
953  }
954 }
955 
956 
957 
959 {
960  // Make sure there is no leftover information
961  this->elem_cmap_elem_ids.clear();
962  this->elem_cmap_side_ids.clear();
963  this->elem_cmap_proc_ids.clear();
964 
965  // Allocate enough space for all our element maps
966  this->elem_cmap_elem_ids.resize(this->num_elem_cmaps);
967  this->elem_cmap_side_ids.resize(this->num_elem_cmaps);
968  this->elem_cmap_proc_ids.resize(this->num_elem_cmaps);
969  {
970  unsigned cnt=0; // Index into vectors
972  it = this->proc_border_elem_sets.begin(),
973  end = this->proc_border_elem_sets.end();
974 
975  for (; it != end; ++it)
976  {
977  // Make sure the current elem_cmap_id matches the index in our map of node intersections
978  libmesh_assert_equal_to (static_cast<unsigned>(this->elem_cmap_ids[cnt]), it->first);
979 
980  // Get reference to the set of IDs to be packed into the vector
981  std::set<std::pair<unsigned,unsigned>> & elem_set = it->second;
982 
983  // Resize the vectors to receive their payload
984  this->elem_cmap_elem_ids[cnt].resize(elem_set.size());
985  this->elem_cmap_side_ids[cnt].resize(elem_set.size());
986  this->elem_cmap_proc_ids[cnt].resize(elem_set.size());
987 
988  std::set<std::pair<unsigned,unsigned>>::iterator elem_set_iter = elem_set.begin();
989 
990  // Pack the vectors with elem IDs, side IDs, and processor IDs.
991  for (std::size_t j=0, eceis=this->elem_cmap_elem_ids[cnt].size(); j<eceis; ++j, ++elem_set_iter)
992  {
993  this->elem_cmap_elem_ids[cnt][j] =
994  libmesh_map_find(libmesh_elem_num_to_exodus, elem_set_iter->first);
995  this->elem_cmap_side_ids[cnt][j] = elem_set_iter->second; // Side ID, this has already been converted above
996  this->elem_cmap_proc_ids[cnt][j] = it->first; // All have the same processor ID
997  }
998 
999  // increment vector index to go to next processor
1000  cnt++;
1001  }
1002  } // end scope for packing
1003 }
1004 
1005 
1006 
1007 
1008 
1010 {
1011  // Make sure we don't have any leftover information
1012  this->node_mapi.clear();
1013  this->node_mapb.clear();
1014  this->node_mape.clear();
1015 
1016  // Make sure there's enough space to hold all our node IDs
1017  this->node_mapi.resize(this->internal_node_ids.size());
1018  this->node_mapb.resize(this->border_node_ids.size());
1019 
1020  // Copy set contents into vectors
1021  {
1022  unsigned cnt = 0;
1023  for (const auto & id : this->internal_node_ids)
1024  this->node_mapi[cnt++] = libmesh_map_find(libmesh_node_num_to_exodus, id);
1025  }
1026 
1027  {
1028  unsigned cnt=0;
1029  for (const auto & id : this->border_node_ids)
1030  this->node_mapb[cnt++] = libmesh_map_find(libmesh_node_num_to_exodus, id);
1031  }
1032 }
1033 
1034 
1035 
1036 
1037 
1039 {
1040  // Make sure there's no left-over information
1041  this->node_cmap_node_ids.clear();
1042  this->node_cmap_proc_ids.clear();
1043 
1044  libmesh_assert_less_equal
1045  (this->proc_nodes_touched_intersections.size(),
1046  std::size_t(this->num_node_cmaps));
1047 
1048  // Allocate enough space for all our node maps
1049  this->node_cmap_node_ids.resize(this->num_node_cmaps);
1050  this->node_cmap_proc_ids.resize(this->num_node_cmaps);
1051  {
1052  unsigned cnt=0; // Index into vectors
1054  it = this->proc_nodes_touched_intersections.begin(),
1055  end = this->proc_nodes_touched_intersections.end();
1056 
1057  for (; it != end; ++it)
1058  {
1059  // Make sure the current node_cmap_id matches the index in our map of node intersections
1060  libmesh_assert_equal_to (static_cast<unsigned>(this->node_cmap_ids[cnt]), it->first);
1061 
1062  // Get reference to the set of IDs to be packed into the vector.
1063  std::set<unsigned> & node_set = it->second;
1064 
1065  // Resize the vectors to receive their payload
1066  this->node_cmap_node_ids[cnt].resize(node_set.size());
1067  this->node_cmap_proc_ids[cnt].resize(node_set.size());
1068 
1069  std::set<unsigned>::iterator node_set_iter = node_set.begin();
1070 
1071  // Pack the vectors with node IDs and processor IDs.
1072  for (std::size_t j=0, nceis=this->node_cmap_node_ids[cnt].size(); j<nceis; ++j, ++node_set_iter)
1073  {
1074  this->node_cmap_node_ids[cnt][j] =
1075  libmesh_map_find(libmesh_node_num_to_exodus, *node_set_iter);
1076  this->node_cmap_proc_ids[cnt][j] = it->first;
1077  }
1078 
1079  // increment vector index to go to next processor
1080  cnt++;
1081  }
1082  } // end scope for packing
1083 
1084  // Print out the vectors we just packed
1085  if (verbose)
1086  {
1087  for (auto i : index_range(this->node_cmap_node_ids))
1088  {
1089  libMesh::out << "[" << this->processor_id() << "] nodes communicated to proc "
1090  << this->node_cmap_ids[i]
1091  << " = ";
1092  for (const auto & node_id : this->node_cmap_node_ids[i])
1093  libMesh::out << node_id << " ";
1094  libMesh::out << std::endl;
1095  }
1096 
1097  for (const auto & id_vec : this->node_cmap_node_ids)
1098  {
1099  libMesh::out << "[" << this->processor_id() << "] processor ID node communicated to = ";
1100  for (const auto & proc_id : id_vec)
1101  libMesh::out << proc_id << " ";
1102  libMesh::out << std::endl;
1103  }
1104  }
1105 }
1106 
1107 
1108 
1109 
1111 {
1112  // For the nodes, these are the number of entries in the sets in proc_nodes_touched_intersections
1113  // map computed above. Note: this map does not contain self-intersections so we can loop over it
1114  // directly.
1115  this->node_cmap_node_cnts.clear(); // Make sure we don't have any leftover information...
1116  this->node_cmap_ids.clear(); // Make sure we don't have any leftover information...
1117  this->node_cmap_node_cnts.resize(this->num_node_cmaps);
1118  this->node_cmap_ids.resize(this->num_node_cmaps);
1119 
1120  {
1121  unsigned cnt=0; // Index into the vector
1123  it = this->proc_nodes_touched_intersections.begin(),
1124  end = this->proc_nodes_touched_intersections.end();
1125 
1126  for (; it != end; ++it)
1127  {
1128  this->node_cmap_ids[cnt] = it->first; // The ID of the proc we communicate with
1129  this->node_cmap_node_cnts[cnt] = cast_int<int>(it->second.size()); // The number of nodes we communicate
1130  cnt++; // increment vector index!
1131  }
1132  }
1133 
1134  // Print the packed vectors we just filled
1135  if (verbose)
1136  {
1137  libMesh::out << "[" << this->processor_id() << "] node_cmap_node_cnts = ";
1138  for (const auto & node_cnt : node_cmap_node_cnts)
1139  libMesh::out << node_cnt << ", ";
1140  libMesh::out << std::endl;
1141 
1142  libMesh::out << "[" << this->processor_id() << "] node_cmap_ids = ";
1143  for (const auto & node_id : node_cmap_ids)
1144  libMesh::out << node_id << ", ";
1145  libMesh::out << std::endl;
1146  }
1147 
1148  // For the elements, we have not yet computed all this information..
1149  this->elem_cmap_elem_cnts.clear(); // Make sure we don't have any leftover information...
1150  this->elem_cmap_ids.clear(); // Make sure we don't have any leftover information...
1151  this->elem_cmap_elem_cnts.resize(this->num_elem_cmaps);
1152  this->elem_cmap_ids.resize(this->num_elem_cmaps);
1153 
1154  // Pack the elem_cmap_ids and elem_cmap_elem_cnts vectors
1155  {
1156  unsigned cnt=0; // Index into the vectors we're filling
1158  it = this->proc_border_elem_sets.begin(),
1159  end = this->proc_border_elem_sets.end();
1160 
1161  for (; it != end; ++it)
1162  {
1163  this->elem_cmap_ids[cnt] = it->first; // The ID of the proc we communicate with
1164  this->elem_cmap_elem_cnts[cnt] = cast_int<int>(it->second.size()); // The number of elems we communicate to/from that proc
1165  cnt++; // increment vector index!
1166  }
1167  }
1168 
1169  // Print the packed vectors we just filled
1170  if (verbose)
1171  {
1172  libMesh::out << "[" << this->processor_id() << "] elem_cmap_elem_cnts = ";
1173  for (const auto & elem_cnt : elem_cmap_elem_cnts)
1174  libMesh::out << elem_cnt << ", ";
1175  libMesh::out << std::endl;
1176 
1177  libMesh::out << "[" << this->processor_id() << "] elem_cmap_ids = ";
1178  for (const auto & elem_id : elem_cmap_ids)
1179  libMesh::out << elem_id << ", ";
1180  libMesh::out << std::endl;
1181  }
1182 }
1183 
1184 
1185 
1186 
1187 void
1189 {
1190  // Set of all local, active element IDs. After we have identified border element
1191  // IDs, the set_difference between this set and the border_elem_ids set will give us
1192  // the set of internal_elem_ids.
1193  std::set<unsigned> all_elem_ids;
1194 
1195  // A set of processor IDs which elements on this processor have as
1196  // neighbors. The size of this set will determine the number of
1197  // element communication maps in Exodus.
1198  std::set<unsigned> neighboring_processor_ids;
1199 
1200  for (const auto & elem : pmesh.active_local_element_ptr_range())
1201  {
1202  // Add this Elem's ID to all_elem_ids, later we will take the difference
1203  // between this set and the set of border_elem_ids, to get the set of
1204  // internal_elem_ids.
1205  all_elem_ids.insert(elem->id());
1206 
1207  // Will be set to true if element is determined to be a border element
1208  bool is_border_elem = false;
1209 
1210  // Construct a conversion object for this Element. This will help us map
1211  // Libmesh numberings into Nemesis numberings for sides.
1212  const auto & conv = get_conversion(elem->type());
1213 
1214  // Add all this element's node IDs to the set of all node IDs.
1215  // The set of internal_node_ids will be the set difference between
1216  // the set of all nodes and the set of border nodes.
1217  //
1218  // In addition, if any node of a local node is listed in the
1219  // border nodes list, then this element goes into the proc_border_elem_sets.
1220  // Note that there is not a 1:1 correspondence between
1221  // border_elem_ids and the entries which go into proc_border_elem_sets.
1222  // The latter is for communication purposes, ie determining which elements
1223  // should be shared between processors.
1224  for (auto node : elem->node_index_range())
1225  this->nodes_attached_to_local_elems.insert(elem->node_id(node));
1226 
1227  // Loop over element's neighbors, see if it has a neighbor which is off-processor
1228  for (auto n : elem->side_index_range())
1229  {
1230  if (elem->neighbor_ptr(n) != nullptr)
1231  {
1232  unsigned neighbor_proc_id = elem->neighbor_ptr(n)->processor_id();
1233 
1234  // If my neighbor has a different processor ID, I must be a border element.
1235  // Also track the neighboring processor ID if it is are different from our processor ID
1236  if (neighbor_proc_id != this->processor_id())
1237  {
1238  is_border_elem = true;
1239  neighboring_processor_ids.insert(neighbor_proc_id);
1240 
1241  // Convert libmesh side(n) of this element into a side ID for Nemesis
1242  unsigned nemesis_side_id = conv.get_inverse_side_map(n);
1243 
1244  if (verbose)
1245  libMesh::out << "[" << this->processor_id() << "] LibMesh side "
1246  << n
1247  << " mapped to (1-based) Exodus side "
1248  << nemesis_side_id
1249  << std::endl;
1250 
1251  // Add this element's ID and the ID of the side which is on the boundary
1252  // to the set of border elements for this processor.
1253  // Note: if the set does not already exist, this creates it.
1254  this->proc_border_elem_sets[ neighbor_proc_id ].emplace(elem->id(), nemesis_side_id);
1255  }
1256  }
1257  } // end for loop over neighbors
1258 
1259  // If we're on a border element, add it to the set
1260  if (is_border_elem)
1261  this->border_elem_ids.insert( elem->id() );
1262 
1263  } // end for loop over active local elements
1264 
1265  // Take the set_difference between all elements and border elements to get internal
1266  // element IDs
1267  std::set_difference(all_elem_ids.begin(), all_elem_ids.end(),
1268  this->border_elem_ids.begin(), this->border_elem_ids.end(),
1269  std::inserter(this->internal_elem_ids, this->internal_elem_ids.end()));
1270 
1271  // Take the set_difference between all nodes and border nodes to get internal nodes
1272  std::set_difference(this->nodes_attached_to_local_elems.begin(), this->nodes_attached_to_local_elems.end(),
1273  this->border_node_ids.begin(), this->border_node_ids.end(),
1274  std::inserter(this->internal_node_ids, this->internal_node_ids.end()));
1275 
1276  if (verbose)
1277  {
1278  libMesh::out << "[" << this->processor_id() << "] neighboring_processor_ids = ";
1279  for (const auto & id : neighboring_processor_ids)
1280  libMesh::out << id << " ";
1281  libMesh::out << std::endl;
1282  }
1283 
1284  // The size of the neighboring_processor_ids set should be the number of element communication maps
1285  this->num_elem_cmaps =
1286  cast_int<int>(neighboring_processor_ids.size());
1287 
1288  if (verbose)
1289  libMesh::out << "[" << this->processor_id() << "] "
1290  << "Number of neighboring processor IDs="
1291  << this->num_elem_cmaps
1292  << std::endl;
1293 
1294  if (verbose)
1295  {
1296  // Print out counts of border elements for each processor
1297  for (const auto & [proc_id, set] : proc_border_elem_sets)
1298  {
1299  libMesh::out << "[" << this->processor_id() << "] "
1300  << "Proc "
1301  << proc_id << " communicates "
1302  << set.size() << " elements." << std::endl;
1303  }
1304  }
1305 
1306  // Store the number of internal and border elements, and the number of internal nodes,
1307  // to be written to the Nemesis file.
1308  this->num_internal_elems =
1309  cast_int<int>(this->internal_elem_ids.size());
1310  this->num_border_elems =
1311  cast_int<int>(this->border_elem_ids.size());
1312  this->num_internal_nodes =
1313  cast_int<int>(this->internal_node_ids.size());
1314 
1315  if (verbose)
1316  {
1317  libMesh::out << "[" << this->processor_id() << "] num_internal_nodes=" << this->num_internal_nodes << std::endl;
1318  libMesh::out << "[" << this->processor_id() << "] num_border_nodes=" << this->num_border_nodes << std::endl;
1319  libMesh::out << "[" << this->processor_id() << "] num_border_elems=" << this->num_border_elems << std::endl;
1320  libMesh::out << "[" << this->processor_id() << "] num_internal_elems=" << this->num_internal_elems << std::endl;
1321  }
1322 }
1323 
1324 
1325 
1327 {
1328  // 1.) Get reference to the set of side boundary IDs
1329  std::set<boundary_id_type> global_side_boundary_ids
1330  (pmesh.get_boundary_info().get_side_boundary_ids().begin(),
1331  pmesh.get_boundary_info().get_side_boundary_ids().end());
1332 
1333  // 2.) Gather boundary side IDs from other processors
1334  this->comm().set_union(global_side_boundary_ids);
1335 
1336  // 3.) Now global_side_boundary_ids actually contains a global list of all side boundary IDs
1337  this->num_side_sets_global =
1338  cast_int<int>(global_side_boundary_ids.size());
1339 
1340  // 4.) Pack these sidesets into a vector so they can be written by Nemesis
1341  this->global_sideset_ids.clear(); // Make sure there is no leftover information
1342  this->global_sideset_ids.insert(this->global_sideset_ids.end(),
1343  global_side_boundary_ids.begin(),
1344  global_side_boundary_ids.end());
1345 
1346  if (verbose)
1347  {
1348  libMesh::out << "[" << this->processor_id() << "] global_sideset_ids = ";
1349  for (const auto & id : this->global_sideset_ids)
1350  libMesh::out << id << ", ";
1351  libMesh::out << std::endl;
1352  }
1353 
1354  // We also need global counts of sides in each of the sidesets.
1355  // Build a list of (elem, side, bc) tuples.
1356  typedef std::tuple<dof_id_type, unsigned short int, boundary_id_type> Tuple;
1357  std::vector<Tuple> bc_triples = pmesh.get_boundary_info().build_side_list();
1358 
1359  // Iterators to the beginning and end of the current range.
1360  std::vector<Tuple>::iterator
1361  it = bc_triples.begin(),
1362  new_end = bc_triples.end();
1363 
1364  while (it != new_end)
1365  {
1366  if (pmesh.elem_ref(std::get<0>(*it)).processor_id() != this->processor_id())
1367  {
1368  // Back up the new end iterators to prepare for swap
1369  --new_end;
1370 
1371  // Swap places, the non-local elem will now be "past-the-end"
1372  std::swap (*it, *new_end);
1373  }
1374  else // elem is local, go to next
1375  ++it;
1376  }
1377 
1378  // Erase from "new" end to old.
1379  bc_triples.erase(new_end, bc_triples.end());
1380 
1381  this->num_global_side_counts.clear(); // Make sure we don't have any leftover information
1382  this->num_global_side_counts.resize(this->global_sideset_ids.size());
1383 
1384  // Get the count for each global sideset ID
1385  for (auto i : index_range(global_sideset_ids))
1386  {
1387  int id = global_sideset_ids[i];
1388  this->num_global_side_counts[i] =
1389  cast_int<int>(std::count_if(bc_triples.begin(),
1390  bc_triples.end(),
1391  [id](const Tuple & t)->bool { return std::get<2>(t) == id; }));
1392  }
1393 
1394  if (verbose)
1395  {
1396  libMesh::out << "[" << this->processor_id() << "] num_global_side_counts = ";
1397  for (const auto & cnt : this->num_global_side_counts)
1398  libMesh::out << cnt << ", ";
1399  libMesh::out << std::endl;
1400  }
1401 
1402  // Finally sum up the result
1403  this->comm().sum(this->num_global_side_counts);
1404 
1405  if (verbose)
1406  {
1407  libMesh::out << "[" << this->processor_id() << "] num_global_side_counts = ";
1408  for (const auto & cnt : this->num_global_side_counts)
1409  libMesh::out << cnt << ", ";
1410  libMesh::out << std::endl;
1411  }
1412 }
1413 
1414 
1415 
1416 
1417 
1418 
1420 {
1421  std::set<boundary_id_type> local_node_boundary_ids;
1422 
1423  // 1.) Get reference to the set of node boundary IDs *for this processor*
1424  std::set<boundary_id_type> global_node_boundary_ids
1425  (pmesh.get_boundary_info().get_node_boundary_ids().begin(),
1426  pmesh.get_boundary_info().get_node_boundary_ids().end());
1427 
1428  // Save a copy of the local_node_boundary_ids...
1429  local_node_boundary_ids = global_node_boundary_ids;
1430 
1431  // 2.) Gather boundary node IDs from other processors
1432  this->comm().set_union(global_node_boundary_ids);
1433 
1434  // 3.) Now global_node_boundary_ids actually contains a global list of all node boundary IDs
1435  this->num_node_sets_global =
1436  cast_int<int>(global_node_boundary_ids.size());
1437 
1438  // 4.) Create a vector<int> from the global_node_boundary_ids set
1439  this->global_nodeset_ids.clear();
1440  this->global_nodeset_ids.insert(this->global_nodeset_ids.end(),
1441  global_node_boundary_ids.begin(),
1442  global_node_boundary_ids.end());
1443 
1444  if (verbose)
1445  {
1446  libMesh::out << "[" << this->processor_id() << "] global_nodeset_ids = ";
1447  for (const auto & id : global_nodeset_ids)
1448  libMesh::out << id << ", ";
1449  libMesh::out << std::endl;
1450 
1451  libMesh::out << "[" << this->processor_id() << "] local_node_boundary_ids = ";
1452  for (const auto & id : local_node_boundary_ids)
1453  libMesh::out << id << ", ";
1454  libMesh::out << std::endl;
1455  }
1456 
1457  // 7.) We also need to know the number of nodes which is in each of the nodesets, globally.
1458 
1459  // Build list of (node-id, bc-id) tuples.
1460  typedef std::tuple<dof_id_type, boundary_id_type> Tuple;
1461  std::vector<Tuple> bc_tuples = pmesh.get_boundary_info().build_node_list();
1462 
1463  if (verbose)
1464  {
1465  libMesh::out << "[" << this->processor_id() << "] boundary_node_list.size()="
1466  << bc_tuples.size() << std::endl;
1467  libMesh::out << "[" << this->processor_id() << "] (boundary_node_id, boundary_id) = ";
1468  for (const auto & t : bc_tuples)
1469  libMesh::out << "(" << std::get<0>(t) << ", " << std::get<1>(t) << ") ";
1470  libMesh::out << std::endl;
1471  }
1472 
1473  // Now get the global information. In this case, we only want to count boundary
1474  // information for nodes *owned* by this processor, so there are no duplicates.
1475 
1476  // Make sure we don't have any left over information
1477  this->num_global_node_counts.clear();
1478  this->num_global_node_counts.resize(this->global_nodeset_ids.size());
1479 
1480  // Unfortunately, we can't just count up all occurrences of a given id,
1481  // that would give us duplicate entries when we do the parallel summation.
1482  // So instead, only count entries for nodes owned by this processor.
1483  // Start by getting rid of all non-local node entries from the vectors.
1484  std::vector<Tuple>::iterator
1485  it = bc_tuples.begin(),
1486  new_end = bc_tuples.end();
1487 
1488  while (it != new_end)
1489  {
1490  if (pmesh.node_ptr(std::get<0>(*it))->processor_id() != this->processor_id())
1491  {
1492  // Back up the new end iterators to prepare for swap
1493  --new_end;
1494 
1495  // Swap places, the non-local node will now be "past-the-end"
1496  std::swap(*it, *new_end);
1497  }
1498  else // node is local, go to next
1499  ++it;
1500  }
1501 
1502  // Erase from "new" end to old end.
1503  bc_tuples.erase(new_end, bc_tuples.end());
1504 
1505  // Now we can do the local count for each ID...
1506  for (auto i : index_range(global_nodeset_ids))
1507  {
1508  int id = this->global_nodeset_ids[i];
1509  this->num_global_node_counts[i] =
1510  cast_int<int>(std::count_if(bc_tuples.begin(),
1511  bc_tuples.end(),
1512  [id](const Tuple & t)->bool { return std::get<1>(t) == id; }));
1513  }
1514 
1515  // And finally we can sum them up
1516  this->comm().sum(this->num_global_node_counts);
1517 
1518  if (verbose)
1519  {
1520  libMesh::out << "[" << this->processor_id() << "] num_global_node_counts = ";
1521  for (const auto & cnt : num_global_node_counts)
1522  libMesh::out << cnt << ", ";
1523  libMesh::out << std::endl;
1524  }
1525 }
1526 
1527 
1528 
1529 
1531 {
1532  // 1.) Loop over active local elements, build up set of subdomain IDs.
1533  std::set<subdomain_id_type> global_subdomain_ids;
1534 
1535  // This map keeps track of the number of elements in each subdomain over all processors
1536  std::map<subdomain_id_type, unsigned> global_subdomain_counts;
1537 
1538  for (const auto & elem : pmesh.active_local_element_ptr_range())
1539  {
1540  subdomain_id_type cur_subdomain = elem->subdomain_id();
1541 
1542  /*
1543  // We can't have a zero subdomain ID in Exodus (for some reason?)
1544  // so map zero subdomains to a max value...
1545  if (cur_subdomain == 0)
1546  cur_subdomain = std::numeric_limits<subdomain_id_type>::max();
1547  */
1548 
1549  global_subdomain_ids.insert(cur_subdomain);
1550 
1551  // Increment the count of elements in this subdomain
1552  global_subdomain_counts[cur_subdomain]++;
1553  }
1554 
1555  // We're next going to this->comm().sum the subdomain counts, so save the local counts
1556  this->local_subdomain_counts = global_subdomain_counts;
1557 
1558  {
1559  // 2.) Copy local subdomain IDs into a vector for communication
1560  std::vector<subdomain_id_type> global_subdomain_ids_vector(global_subdomain_ids.begin(),
1561  global_subdomain_ids.end());
1562 
1563  // 3.) Gather them into an enlarged vector
1564  this->comm().allgather(global_subdomain_ids_vector);
1565 
1566  // 4.) Insert any new IDs into the set (any duplicates will be dropped)
1567  global_subdomain_ids.insert(global_subdomain_ids_vector.begin(),
1568  global_subdomain_ids_vector.end());
1569  }
1570 
1571  // 5.) Now global_subdomain_ids actually contains a global list of all subdomain IDs
1572  this->num_elem_blks_global =
1573  cast_int<int>(global_subdomain_ids.size());
1574 
1575  // Print the number of elements found locally in each subdomain
1576  if (verbose)
1577  {
1578  libMesh::out << "[" << this->processor_id() << "] ";
1579  for (const auto & [subdomain_id, cnt] : global_subdomain_counts)
1580  {
1581  libMesh::out << "ID: "
1582  << static_cast<unsigned>(subdomain_id)
1583  << ", Count: " << cnt << ", ";
1584  }
1585  libMesh::out << std::endl;
1586  }
1587 
1588  // 6.) this->comm().sum up the number of elements in each block. We know the global
1589  // subdomain IDs, so pack them into a vector one by one. Use a vector of int since
1590  // that is what Nemesis wants
1591  this->global_elem_blk_cnts.resize(global_subdomain_ids.size());
1592 
1593  unsigned cnt=0;
1594  // Find the entry in the local map, note: if not found, will be created with 0 default value, which is OK...
1595  for (const auto & id : global_subdomain_ids)
1596  this->global_elem_blk_cnts[cnt++] = global_subdomain_counts[id];
1597 
1598  // Sum up subdomain counts from all processors
1599  this->comm().sum(this->global_elem_blk_cnts);
1600 
1601  if (verbose)
1602  {
1603  libMesh::out << "[" << this->processor_id() << "] global_elem_blk_cnts = ";
1604  for (const auto & bc : this->global_elem_blk_cnts)
1605  libMesh::out << bc << ", ";
1606  libMesh::out << std::endl;
1607  }
1608 
1609  // 7.) Create a vector<int> from the global_subdomain_ids set, for passing to Nemesis
1610  this->global_elem_blk_ids.clear();
1611  this->global_elem_blk_ids.insert(this->global_elem_blk_ids.end(), // pos
1612  global_subdomain_ids.begin(),
1613  global_subdomain_ids.end());
1614 
1615  if (verbose)
1616  {
1617  libMesh::out << "[" << this->processor_id() << "] global_elem_blk_ids = ";
1618  for (const auto & id : this->global_elem_blk_ids)
1619  libMesh::out << id << ", ";
1620  libMesh::out << std::endl;
1621  }
1622 
1623 
1624  // 8.) We will call put_eb_info_global later, it must be called after this->put_init_global().
1625 }
1626 
1627 
1628 
1629 
1631 {
1632  // If we don't have any local subdomains, it had better be because
1633  // we don't have any local elements
1634 #ifdef DEBUG
1635  if (local_subdomain_counts.empty())
1636  {
1637  libmesh_assert(pmesh.active_local_elements_begin() ==
1638  pmesh.active_local_elements_end());
1640  }
1641 #endif
1642 
1643  // Elements have to be numbered contiguously based on what block
1644  // number they are in. Therefore we have to do a bit of work to get
1645  // the block (ie subdomain) numbers first and store them off as
1646  // block_ids.
1647 
1648  // Make sure there is no leftover information in the subdomain_map, and reserve
1649  // enough space to store the elements we need.
1650  this->subdomain_map.clear();
1651  for (const auto & [sbd_id, cnt] : local_subdomain_counts)
1652  {
1653  if (verbose)
1654  {
1655  libMesh::out << "[" << this->processor_id() << "] "
1656  << "local_subdomain_counts [" << static_cast<unsigned>(sbd_id) << "]= "
1657  << cnt
1658  << std::endl;
1659  }
1660 
1661  this->subdomain_map[sbd_id].reserve(cnt);
1662  }
1663 
1664 
1665  // First loop over the elements to figure out which elements are in which subdomain
1666  for (const auto & elem : pmesh.active_local_element_ptr_range())
1667  {
1668  // Grab the nodes while we're here.
1669  for (auto n : elem->node_index_range())
1670  this->nodes_attached_to_local_elems.insert( elem->node_id(n) );
1671 
1672  subdomain_id_type cur_subdomain = elem->subdomain_id();
1673 
1674  this->subdomain_map[cur_subdomain].push_back(elem->id());
1675  }
1676 
1677  // Set num_nodes which is used by exodusII_io_helper
1678  this->num_nodes =
1679  cast_int<int>(this->nodes_attached_to_local_elems.size());
1680 
1681  // Now come up with a 1-based numbering for these nodes
1682  this->exodus_node_num_to_libmesh.clear(); // Make sure it's empty
1683  this->exodus_node_num_to_libmesh.reserve(this->nodes_attached_to_local_elems.size());
1684 
1685  // Also make sure there's no leftover information in the map which goes the
1686  // other direction.
1687  this->libmesh_node_num_to_exodus.clear();
1688 
1689  // Set the map for nodes
1690  for (const auto & id : nodes_attached_to_local_elems)
1691  {
1692  // I.e. given exodus_node_id,
1693  // exodus_node_num_to_libmesh[ exodus_node_id ] returns the
1694  // libmesh ID for that node, plus one.
1695  // Here we index libMesh IDs with an offset of 1 because they're
1696  // the literal numbers that get written to the exodus file, but
1697  // we index Exodus IDs with an offset of 0 because we read that
1698  // exodus data into a C++ vector. Confused yet?
1699  this->exodus_node_num_to_libmesh.push_back(id+1);
1700 
1701  // Likewise, given libmesh_node_id,
1702  // libmesh_node_num_to_exodus[ libmesh_node_id] returns the
1703  // *Exodus* ID for that node. Unlike the
1704  // exodus_node_num_to_libmesh vector above, this one is a
1705  // std::map. We're never handing a data buffer from it over to
1706  // another API so we don't need to do any weird offsets with it.
1707  this->libmesh_node_num_to_exodus[id] =
1708  this->exodus_node_num_to_libmesh.size(); // should never be zero...
1709  }
1710 
1711  // Now we're going to loop over the subdomain map and build a few things right
1712  // now that we'll use later.
1713 
1714  // First make sure our data structures don't have any leftover data...
1715  this->exodus_elem_num_to_libmesh.clear();
1716  this->block_ids.clear();
1717  this->libmesh_elem_num_to_exodus.clear();
1718 
1719  // Now loop over each subdomain and get a unique numbering for the elements
1720  for (auto & [block_id, elem_ids_this_subdomain] : subdomain_map)
1721  {
1722  block_ids.push_back(block_id);
1723 
1724  // The code below assumes this subdomain block is not empty, make sure that's the case!
1725  libmesh_error_msg_if(elem_ids_this_subdomain.size() == 0,
1726  "Error, no element IDs found in subdomain " << block_id);
1727 
1728  // Use the first element in this block to get representative information.
1729  // Note that Exodus assumes all elements in a block are of the same type!
1730  // We are using that same assumption here!
1731  const auto & conv = get_conversion
1732  (pmesh.elem_ref(elem_ids_this_subdomain[0]).type());
1733  this->num_nodes_per_elem =
1734  pmesh.elem_ref(elem_ids_this_subdomain[0]).n_nodes();
1735 
1736  // Get a reference to the connectivity vector for this subdomain. This vector
1737  // is most likely empty, we are going to fill it up now.
1738  std::vector<int> & current_block_connectivity = this->block_id_to_elem_connectivity[block_id];
1739 
1740  // Just in case it's not already empty...
1741  current_block_connectivity.clear();
1742  current_block_connectivity.resize(elem_ids_this_subdomain.size() * this->num_nodes_per_elem);
1743 
1744  for (auto i : index_range(elem_ids_this_subdomain))
1745  {
1746  auto elem_id = elem_ids_this_subdomain[i];
1747 
1748  // Set the number map for elements
1749  // exodus_elem_num_to_libmesh[ exodus_node_id ] returns the
1750  // libmesh ID for that element, plus one.
1751  // Like with nodes above, we index libMesh IDs with an
1752  // offset of 1 because they're the literal numbers that get
1753  // written to the exodus file, but we index Exodus IDs with
1754  // an offset of 0 because we read that exodus data into a
1755  // C++ vector.
1756  this->exodus_elem_num_to_libmesh.push_back(elem_id+1);
1757 
1758  // Likewise, given libmesh elem_id,
1759  // libmesh_elem_num_to_exodus[ elem_id ] returns the
1760  // *Exodus* ID for that node. Unlike the
1761  // exodus_elem_num_to_libmesh vector above, this one is a
1762  // std::map. We're never handing a data buffer from it over to
1763  // another API so we don't need to do any weird offsets with it.
1764  this->libmesh_elem_num_to_exodus[elem_id] =
1765  this->exodus_elem_num_to_libmesh.size();
1766 
1767  const Elem & elem = pmesh.elem_ref(elem_id);
1768 
1769  // Exodus/Nemesis want every block to have the same element type
1770  // libmesh_assert_equal_to (elem->type(), conv.libmesh_elem_type());
1771 
1772  // But we can get away with writing e.g. HEX8 and INFHEX8 in
1773  // the same block...
1774  libmesh_assert_equal_to (elem.n_nodes(), Elem::build(conv.libmesh_elem_type(), nullptr)->n_nodes());
1775 
1776  for (auto j : make_range(this->num_nodes_per_elem))
1777  {
1778  const unsigned int connect_index = (i*this->num_nodes_per_elem)+j;
1779  const unsigned int elem_node_index = conv.get_inverse_node_map(j); // inverse node map is used for writing
1780 
1781  current_block_connectivity[connect_index] =
1782  libmesh_map_find(libmesh_node_num_to_exodus,
1783  elem.node_id(elem_node_index));
1784  }
1785  } // End loop over elems in this subdomain
1786  } // end loop over subdomain_map
1787 }
1788 
1789 
1790 
1791 
1792 
1794 {
1795  // The set which will eventually contain the IDs of "border nodes". These are nodes
1796  // that lie on the boundary between one or more processors.
1797  //std::set<unsigned> border_node_ids;
1798 
1799  // map from processor ID to set of nodes which elements from this processor "touch",
1800  // that is,
1801  // proc_nodes_touched[p] = (set all node IDs found in elements owned by processor p)
1802  std::map<unsigned, std::set<unsigned>> proc_nodes_touched;
1803 
1804 
1805  // We are going to create a lot of intermediate data structures here, so make sure
1806  // as many as possible all cleaned up by creating scope!
1807  {
1808  // Loop over active (not just active local) elements, make sets of node IDs for each
1809  // processor which has an element that "touches" a node.
1810  for (const auto & elem : pmesh.active_element_ptr_range())
1811  {
1812  // Get reference to the set for this processor. If it does not exist
1813  // it will be created.
1814  std::set<unsigned> & set_p = proc_nodes_touched[ elem->processor_id() ];
1815 
1816  // Insert all nodes touched by this element into the set
1817  for (auto node : elem->node_index_range())
1818  set_p.insert(elem->node_id(node));
1819  }
1820 
1821  if (verbose)
1822  {
1823  libMesh::out << "[" << this->processor_id()
1824  << "] proc_nodes_touched contains "
1825  << proc_nodes_touched.size()
1826  << " sets of nodes."
1827  << std::endl;
1828 
1829  for (const auto & [proc_id, set] : proc_nodes_touched)
1830  libMesh::out << "[" << this->processor_id()
1831  << "] proc_nodes_touched[" << proc_id << "] has "
1832  << set.size()
1833  << " entries."
1834  << std::endl;
1835  }
1836 
1837 
1838  // Loop over all the sets we just created and compute intersections with the
1839  // this processor's set. Obviously, don't intersect with ourself.
1840  this->proc_nodes_touched_intersections.clear();
1841  for (auto & [proc_id, other_set] : proc_nodes_touched)
1842  {
1843  // Don't compute intersections with ourself
1844  if (proc_id == this->processor_id())
1845  continue;
1846 
1847  std::set<unsigned int> this_intersection;
1848 
1849  // Otherwise, compute intersection with other processor and ourself
1850  std::set<unsigned> & my_set = proc_nodes_touched[this->processor_id()];
1851 
1852  std::set_intersection(my_set.begin(), my_set.end(),
1853  other_set.begin(), other_set.end(),
1854  std::inserter(this_intersection, this_intersection.end()));
1855 
1856  if (!this_intersection.empty())
1857  this->proc_nodes_touched_intersections.emplace
1858  (proc_id, std::move(this_intersection));
1859  }
1860 
1861  if (verbose)
1862  {
1863  for (const auto & [proc_id, set] : proc_nodes_touched_intersections)
1864  libMesh::out << "[" << this->processor_id()
1865  << "] this->proc_nodes_touched_intersections[" << proc_id << "] has "
1866  << set.size()
1867  << " entries."
1868  << std::endl;
1869  }
1870 
1871  // The number of node communication maps is the number of other processors
1872  // with which we share nodes.
1873  this->num_node_cmaps =
1874  cast_int<int>(proc_nodes_touched_intersections.size());
1875 
1876  // We can't be connecting to more processors than exist outside
1877  // ourselves
1878  libmesh_assert_less (this->num_node_cmaps, this->n_processors());
1879 
1880  // Compute the set_union of all the preceding intersections. This will be the set of
1881  // border node IDs for this processor.
1882  for (auto & pr : proc_nodes_touched_intersections)
1883  {
1884  std::set<unsigned> & other_set = pr.second;
1885  std::set<unsigned> intermediate_result; // Don't think we can insert into one of the sets we're unioning...
1886 
1887  std::set_union(this->border_node_ids.begin(), this->border_node_ids.end(),
1888  other_set.begin(), other_set.end(),
1889  std::inserter(intermediate_result, intermediate_result.end()));
1890 
1891  // Swap our intermediate result into the final set
1892  this->border_node_ids.swap(intermediate_result);
1893  }
1894 
1895  libmesh_assert_less_equal
1896  (this->proc_nodes_touched_intersections.size(),
1897  std::size_t(this->num_node_cmaps));
1898 
1899  if (verbose)
1900  {
1901  libMesh::out << "[" << this->processor_id()
1902  << "] border_node_ids.size()=" << this->border_node_ids.size()
1903  << std::endl;
1904  }
1905  } // end scope for border node ID creation
1906 
1907  // Store the number of border node IDs to be written to Nemesis file
1908  this->num_border_nodes = cast_int<int>(this->border_node_ids.size());
1909 }
1910 
1911 
1912 
1913 
1914 
1916 {
1917  // Write the nodesets. In Nemesis, the idea is to "create space" for the global
1918  // set of boundary nodesets, but to only write node IDs which are local to the current
1919  // processor. This is what is done in Nemesis files created by the "loadbal" script.
1920 
1921  // Store a map of vectors for boundary node IDs on this processor.
1922  // Use a vector of int here so it can be passed directly to Exodus.
1923  std::map<boundary_id_type, std::vector<int>> local_node_boundary_id_lists;
1924 
1925  // FIXME: We should build this list only one time!! We already built it above, but we
1926  // did not have the libmesh to exodus node mapping at that time... for now we'll just
1927  // build it here again, hopefully it's small relative to the size of the entire mesh.
1928 
1929  // Build list of (node-id, bc-id) tuples.
1930  typedef std::tuple<dof_id_type, boundary_id_type> Tuple;
1931  std::vector<Tuple> bc_tuples = mesh.get_boundary_info().build_node_list();
1932 
1933  if (verbose)
1934  {
1935  libMesh::out << "[" << this->processor_id() << "] boundary_node_list.size()="
1936  << bc_tuples.size() << std::endl;
1937  libMesh::out << "[" << this->processor_id() << "] (boundary_node_id, boundary_id) = ";
1938  for (const auto & t : bc_tuples)
1939  libMesh::out << "(" << std::get<0>(t) << ", " << std::get<1>(t) << ") ";
1940  libMesh::out << std::endl;
1941  }
1942 
1943  // For each node in the node list, add it to the vector of node IDs for that
1944  // set for the local processor. This will be used later when writing Exodus
1945  // nodesets.
1946  for (const auto & t : bc_tuples)
1947  {
1948  // Don't try to grab a reference to the vector unless the current node is attached
1949  // to a local element. Otherwise, another processor will be responsible for writing it in its nodeset.
1950  if (const auto it = this->libmesh_node_num_to_exodus.find(std::get<0>(t));
1951  it != this->libmesh_node_num_to_exodus.end())
1952  {
1953  // Get reference to the vector where this node ID will be inserted. If it
1954  // doesn't yet exist, this will create it.
1955  std::vector<int> & current_id_set = local_node_boundary_id_lists[std::get<1>(t)];
1956 
1957  // Push back Exodus-mapped node ID for this set
1958  // TODO: reserve space in these vectors somehow.
1959  current_id_set.push_back( it->second );
1960  }
1961  }
1962 
1963  // See what we got
1964  if (verbose)
1965  {
1966  for (const auto & [bndry_id, set] : local_node_boundary_id_lists)
1967  {
1968  libMesh::out << "[" << this->processor_id() << "] ID: " << bndry_id << ", ";
1969 
1970  // Libmesh node ID (Exodus Node ID)
1971  for (const auto & id : set)
1972  libMesh::out << id << ", ";
1973  libMesh::out << std::endl;
1974  }
1975  }
1976 
1977  // Loop over *global* nodeset IDs, call the Exodus API. Note that some nodesets may be empty
1978  // for a given processor.
1979  if (global_nodeset_ids.size() > 0)
1980  {
1981  NamesData names_table(global_nodeset_ids.size(), MAX_STR_LENGTH);
1982 
1983  for (const auto & nodeset_id : this->global_nodeset_ids)
1984  {
1985  const std::string & current_ns_name =
1986  mesh.get_boundary_info().get_nodeset_name
1987  (cast_int<boundary_id_type>(nodeset_id));
1988 
1989  // Store this name in a data structure that will be used to
1990  // write sideset names to file.
1991  names_table.push_back_entry(current_ns_name);
1992 
1993  if (verbose)
1994  {
1995  libMesh::out << "[" << this->processor_id()
1996  << "] Writing out Exodus nodeset info for ID: " << nodeset_id
1997  << ", Name: " << current_ns_name
1998  << std::endl;
1999  }
2000 
2001  // Convert current global_nodeset_id into an exodus ID, which can't be zero...
2002  int exodus_id = nodeset_id;
2003 
2004  /*
2005  // Exodus can't handle zero nodeset IDs (?) Use max short here since
2006  // when libmesh reads it back in, it will want to store it as a short...
2007  if (exodus_id==0)
2008  exodus_id = std::numeric_limits<short>::max();
2009  */
2010 
2011  // Try to find this boundary ID in the local list we created
2012  if (const auto it = local_node_boundary_id_lists.find (cast_int<boundary_id_type>(nodeset_id));
2013  it == local_node_boundary_id_lists.end())
2014  {
2015  // No nodes found for this boundary ID on this processor
2016  if (verbose)
2017  libMesh::out << "[" << this->processor_id()
2018  << "] No nodeset data for ID: " << nodeset_id
2019  << " on this processor." << std::endl;
2020 
2021  // Call the Exodus interface to write the parameters of this node set
2022  this->ex_err = exII::ex_put_node_set_param(this->ex_id,
2023  exodus_id,
2024  0, /* No nodes for this ID */
2025  0 /* No distribution factors */);
2026  EX_CHECK_ERR(this->ex_err, "Error writing nodeset parameters in Nemesis");
2027 
2028  }
2029  else // Boundary ID *was* found in list
2030  {
2031  // Get reference to the vector of node IDs
2032  const std::vector<int> & current_nodeset_ids = it->second;
2033 
2034  // Call the Exodus interface to write the parameters of this node set
2035  this->ex_err = exII::ex_put_node_set_param(this->ex_id,
2036  exodus_id,
2037  current_nodeset_ids.size(),
2038  0 /* No distribution factors */);
2039 
2040  EX_CHECK_ERR(this->ex_err, "Error writing nodeset parameters in Nemesis");
2041 
2042  // Call Exodus interface to write the actual node IDs for this boundary ID
2043  this->ex_err = exII::ex_put_node_set(this->ex_id,
2044  exodus_id,
2045  current_nodeset_ids.data());
2046 
2047  EX_CHECK_ERR(this->ex_err, "Error writing nodesets in Nemesis");
2048 
2049  }
2050  } // end loop over global nodeset IDs
2051 
2052  // Write out the nodeset names
2053  ex_err = exII::ex_put_names(ex_id,
2054  exII::EX_NODE_SET,
2055  names_table.get_char_star_star());
2056  EX_CHECK_ERR(ex_err, "Error writing nodeset names");
2057  } // end for loop over global nodeset IDs
2058 }
2059 
2060 
2061 
2062 
2064 {
2065  // Write the sidesets. In Nemesis, the idea is to "create space" for the global
2066  // set of boundary sidesets, but to only write sideset IDs which are local to the current
2067  // processor. This is what is done in Nemesis files created by the "loadbal" script.
2068  // See also: ExodusII_IO_Helper::write_sidesets()...
2069 
2070 
2071  // Store a map of vectors for boundary side IDs on this processor.
2072  // Use a vector of int here so it can be passed directly to Exodus.
2073  std::map<boundary_id_type, std::vector<int>> local_elem_boundary_id_lists;
2074  std::map<boundary_id_type, std::vector<int>> local_elem_boundary_id_side_lists;
2075 
2076  // FIXME: We already built this list once, we should reuse that information!
2077  std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> bndry_elem_side_id_list =
2078  mesh.get_boundary_info().build_side_list();
2079 
2080  // Integer looping, skipping non-local elements
2081  for (const auto & t : bndry_elem_side_id_list)
2082  {
2083  // Get pointer to current Elem
2084  const Elem * elem = mesh.elem_ptr(std::get<0>(t));
2085 
2086  std::vector<const Elem *> family;
2087 #ifdef LIBMESH_ENABLE_AMR
2088  // We need to build up active elements if AMR is enabled and add
2089  // them to the exodus sidesets instead of the potentially inactive "parent" elements
2090  // Technically we don't need to "reset" the tree since the vector was just created.
2091  elem->active_family_tree_by_side(family, std::get<1>(t), /*reset tree=*/false);
2092 #else
2093  // If AMR is not even enabled, just push back the element itself
2094  family.push_back( elem );
2095 #endif
2096 
2097  // Loop over all the elements in the family tree, store their converted IDs
2098  // and side IDs to the map's vectors. TODO: Somehow reserve enough space for these
2099  // push_back's...
2100  for (const auto & tree_elem : family)
2101  {
2102  const dof_id_type f_id = tree_elem->id();
2103  const Elem & f = mesh.elem_ref(f_id);
2104 
2105  // If element is local, process it
2106  if (f.processor_id() == this->processor_id())
2107  {
2108  const auto & conv = get_conversion(f.type());
2109 
2110  // Use the libmesh to exodus data structure map to get the proper sideset IDs
2111  // The data structure contains the "collapsed" contiguous ids.
2112  //
2113  // We know the parent element is local, but let's be absolutely sure that all the children have been
2114  // actually mapped to Exodus IDs before we blindly try to add them...
2115  local_elem_boundary_id_lists[ std::get<2>(t) ].push_back( libmesh_map_find(libmesh_elem_num_to_exodus, f_id) );
2116  local_elem_boundary_id_side_lists[ std::get<2>(t) ].push_back(conv.get_inverse_side_map( std::get<1>(t) ));
2117  }
2118  }
2119  }
2120 
2121 
2122  // Loop over *global* sideset IDs, call the Exodus API. Note that some sidesets may be empty
2123  // for a given processor.
2124  if (global_sideset_ids.size() > 0)
2125  {
2126  NamesData names_table(global_sideset_ids.size(), MAX_STR_LENGTH);
2127 
2128  for (const auto & exodus_id : this->global_sideset_ids)
2129  {
2130  const std::string & current_ss_name =
2131  mesh.get_boundary_info().get_sideset_name
2132  (cast_int<boundary_id_type>(exodus_id));
2133 
2134  // Store this name in a data structure that will be used to
2135  // write sideset names to file.
2136  names_table.push_back_entry(current_ss_name);
2137 
2138  if (verbose)
2139  {
2140  libMesh::out << "[" << this->processor_id()
2141  << "] Writing out Exodus sideset info for ID: " << exodus_id
2142  << ", Name: " << current_ss_name
2143  << std::endl;
2144  }
2145 
2146  // Try to find this boundary ID in the local list we created
2147  if (const auto it = local_elem_boundary_id_lists.find (cast_int<boundary_id_type>(exodus_id));
2148  it == local_elem_boundary_id_lists.end())
2149  {
2150  // No sides found for this boundary ID on this processor
2151  if (verbose)
2152  libMesh::out << "[" << this->processor_id()
2153  << "] No sideset data for ID: " << exodus_id
2154  << " on this processor." << std::endl;
2155 
2156  // Call the Exodus interface to write the parameters of this side set
2157  this->ex_err = exII::ex_put_side_set_param(this->ex_id,
2158  exodus_id,
2159  0, /* No sides for this ID */
2160  0 /* No distribution factors */);
2161  EX_CHECK_ERR(this->ex_err, "Error writing sideset parameters in Nemesis");
2162 
2163  }
2164  else // Boundary ID *was* found in list
2165  {
2166  // Get reference to the vector of elem IDs
2167  const std::vector<int> & current_sideset_elem_ids = it->second;
2168 
2169  // Get reference to the vector of side IDs
2170  std::vector<int> & current_sideset_side_ids =
2171  libmesh_map_find(local_elem_boundary_id_side_lists,
2172  cast_int<boundary_id_type>(exodus_id));
2173 
2174  // Call the Exodus interface to write the parameters of this side set
2175  this->ex_err = exII::ex_put_side_set_param(this->ex_id,
2176  exodus_id,
2177  current_sideset_elem_ids.size(),
2178  0 /* No distribution factors */);
2179 
2180  EX_CHECK_ERR(this->ex_err, "Error writing sideset parameters in Nemesis");
2181 
2182  // Call Exodus interface to write the actual side IDs for this boundary ID
2183  this->ex_err = exII::ex_put_side_set(this->ex_id,
2184  exodus_id,
2185  current_sideset_elem_ids.data(),
2186  current_sideset_side_ids.data());
2187 
2188  EX_CHECK_ERR(this->ex_err, "Error writing sidesets in Nemesis");
2189  }
2190  } // end for loop over global sideset IDs
2191 
2192  // Write sideset names to file. Some of these may be blank strings
2193  // if the current processor didn't have all the sideset names for
2194  // any reason...
2195  ex_err = exII::ex_put_names(this->ex_id,
2196  exII::EX_SIDE_SET,
2197  names_table.get_char_star_star());
2198  EX_CHECK_ERR(ex_err, "Error writing sideset names");
2199 
2200  } // end if (global_sideset_ids.size() > 0)
2201 }
2202 
2203 
2204 
2205 void Nemesis_IO_Helper::write_nodal_coordinates(const MeshBase & mesh, bool /*use_discontinuous*/)
2206 {
2207  auto local_num_nodes = this->exodus_node_num_to_libmesh.size();
2208 
2209  x.resize(local_num_nodes);
2210  y.resize(local_num_nodes);
2211  z.resize(local_num_nodes);
2212 
2213  // Just loop over our list outputting the nodes the way we built the map
2214  for (auto i : make_range(local_num_nodes))
2215  {
2216  const Point & pt = mesh.point(this->exodus_node_num_to_libmesh[i]-1);
2217  x[i]=pt(0);
2218  y[i]=pt(1);
2219  z[i]=pt(2);
2220  }
2221 
2222  if (local_num_nodes)
2223  {
2224  // Call Exodus API to write nodal coordinates...
2225  ex_err = exII::ex_put_coord
2226  (ex_id,
2227  x.empty() ? nullptr : MappedOutputVector(x, _single_precision).data(),
2228  y.empty() ? nullptr : MappedOutputVector(y, _single_precision).data(),
2229  z.empty() ? nullptr : MappedOutputVector(z, _single_precision).data());
2230  EX_CHECK_ERR(ex_err, "Error writing node coordinates");
2231 
2232  // And write the nodal map we created for them
2233  ex_err = exII::ex_put_node_num_map(ex_id, this->exodus_node_num_to_libmesh.data());
2234  EX_CHECK_ERR(ex_err, "Error writing node num map");
2235  }
2236  else // Does the Exodus API want us to write empty nodal coordinates?
2237  {
2238  ex_err = exII::ex_put_coord(ex_id, nullptr, nullptr, nullptr);
2239  EX_CHECK_ERR(ex_err, "Error writing empty node coordinates");
2240 
2241  ex_err = exII::ex_put_node_num_map(ex_id, nullptr);
2242  EX_CHECK_ERR(ex_err, "Error writing empty node num map");
2243  }
2244 }
2245 
2246 
2247 
2248 
2249 
2250 void Nemesis_IO_Helper::write_elements(const MeshBase & mesh, bool /*use_discontinuous*/)
2251 {
2252  // Only write elements if there are elements blocks available.
2253  if (this->num_elem_blks_global > 0)
2254  {
2255  // Data structure to store element block names that will be used to
2256  // write the element block names to file.
2257  NamesData names_table(this->num_elem_blks_global, MAX_STR_LENGTH);
2258 
2259  // Loop over all blocks, even if we don't have elements in each block.
2260  // If we don't have elements we need to write out a 0 for that block...
2261  for (auto i : make_range(this->num_elem_blks_global))
2262  {
2263  // Even if there are no elements for this block on the current
2264  // processor, we still want to write its name to file, if
2265  // possible. MeshBase::subdomain_name() will just return an
2266  // empty string if there is no name associated with the current
2267  // block.
2268  names_table.push_back_entry
2269  (mesh.subdomain_name(cast_int<subdomain_id_type>(this->global_elem_blk_ids[i])));
2270 
2271  // Search for the current global block ID in the map
2272  if (const auto it = this->block_id_to_elem_connectivity.find( this->global_elem_blk_ids[i] );
2273  it == this->block_id_to_elem_connectivity.end())
2274  {
2275  // If not found, write a zero to file....
2276  this->ex_err = exII::ex_put_elem_block(this->ex_id,
2277  this->global_elem_blk_ids[i],
2278  "Empty",
2279  0, /* n. elements in this block */
2280  0, /* n. nodes per element */
2281  0); /* number of attributes per element */
2282 
2283  EX_CHECK_ERR(this->ex_err, "Error writing element block from Nemesis.");
2284  }
2285 
2286  // Otherwise, write the actual block information and connectivity to file
2287  else
2288  {
2289  subdomain_id_type block =
2290  cast_int<subdomain_id_type>(it->first);
2291  const std::vector<int> & this_block_connectivity = it->second;
2292  std::vector<dof_id_type> & elements_in_this_block = subdomain_map[block];
2293 
2294  // Use the first element in this block to get representative information.
2295  // Note that Exodus assumes all elements in a block are of the same type!
2296  // We are using that same assumption here!
2297  const auto & conv =
2298  get_conversion(mesh.elem_ref(elements_in_this_block[0]).type());
2299 
2300  this->num_nodes_per_elem =
2301  mesh.elem_ref(elements_in_this_block[0]).n_nodes();
2302 
2303  ex_err = exII::ex_put_elem_block(ex_id,
2304  block,
2305  conv.exodus_elem_type().c_str(),
2306  elements_in_this_block.size(),
2308  0);
2309  EX_CHECK_ERR(ex_err, "Error writing element block from Nemesis.");
2310 
2311  ex_err = exII::ex_put_elem_conn(ex_id,
2312  block,
2313  this_block_connectivity.data());
2314  EX_CHECK_ERR(ex_err, "Error writing element connectivities from Nemesis.");
2315  }
2316  } // end loop over global block IDs
2317 
2318  // Only call this once, not in the loop above!
2319  ex_err = exII::ex_put_elem_num_map(ex_id,
2320  exodus_elem_num_to_libmesh.empty() ? nullptr : exodus_elem_num_to_libmesh.data());
2321  EX_CHECK_ERR(ex_err, "Error writing element map");
2322 
2323  // Write the element block names to file.
2324  ex_err = exII::ex_put_names(ex_id, exII::EX_ELEM_BLOCK, names_table.get_char_star_star());
2325  EX_CHECK_ERR(ex_err, "Error writing element block names");
2326  } // end if (this->num_elem_blks_global > 0)
2327 }
2328 
2329 
2330 
2331 
2332 
2333 void Nemesis_IO_Helper::write_nodal_solution(const std::vector<Number> & values,
2334  const std::vector<std::string> & names,
2335  int timestep)
2336 {
2337  int num_vars = cast_int<int>(names.size());
2338  //int num_values = values.size(); // Not used?
2339 
2340  for (int c=0; c<num_vars; c++)
2341  {
2342 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2343  std::vector<Real> real_parts(num_nodes);
2344  std::vector<Real> imag_parts(num_nodes);
2345  std::vector<Real> magnitudes(num_nodes);
2346 
2347  for (int i=0; i<num_nodes; ++i)
2348  {
2349  Number value = values[(this->exodus_node_num_to_libmesh[i]-1)*num_vars + c];
2350  real_parts[i] = value.real();
2351  imag_parts[i] = value.imag();
2352  magnitudes[i] = std::abs(value);
2353  }
2354  write_nodal_values(3*c+1,real_parts,timestep);
2355  write_nodal_values(3*c+2,imag_parts,timestep);
2356  write_nodal_values(3*c+3,magnitudes,timestep);
2357 #else
2358  std::vector<Number> cur_soln(this->num_nodes);
2359 
2360  // Copy out this variable's solution
2361  for (int i=0; i<this->num_nodes; i++)
2362  cur_soln[i] = values[(this->exodus_node_num_to_libmesh[i]-1)*num_vars + c];
2363 
2364  write_nodal_values(c+1,cur_soln,timestep);
2365 #endif
2366  }
2367 }
2368 
2369 
2370 
2372  const std::vector<std::string> & names,
2373  int timestep,
2374  const std::vector<std::string> & output_names)
2375 {
2376  int num_vars = cast_int<int>(names.size());
2377 
2378  for (int c=0; c<num_vars; c++)
2379  {
2380  // Find the position of names[c] in the output_names vector, if it exists.
2381  auto pos = std::find(output_names.begin(), output_names.end(), names[c]);
2382 
2383  // Skip names[c] if it's not supposed to be output.
2384  if (pos == output_names.end())
2385  continue;
2386 
2387  // Compute the (zero-based) index which determines which
2388  // variable this will be as far as Nemesis is concerned. This
2389  // will be used below in the write_nodal_values() call.
2390  int variable_name_position =
2391  cast_int<int>(std::distance(output_names.begin(), pos));
2392 
2393  // Fill up a std::vector with the dofs for the current variable
2394  std::vector<numeric_index_type> required_indices(this->num_nodes);
2395 
2396  for (int i=0; i<this->num_nodes; i++)
2397  required_indices[i] = static_cast<dof_id_type>(this->exodus_node_num_to_libmesh[i]-1) * num_vars + c;
2398 
2399  // Get the dof values required to write just our local part of
2400  // the solution vector.
2401  std::vector<Number> local_soln;
2402  parallel_soln.localize(local_soln, required_indices);
2403 
2404 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
2405  // Call the ExodusII_IO_Helper function to write the data.
2406  write_nodal_values(variable_name_position + 1, local_soln, timestep);
2407 #else
2408  // We have the local (complex) values. Now extract the real,
2409  // imaginary, and magnitude values from them.
2410  std::vector<Real> real_parts(num_nodes);
2411  std::vector<Real> imag_parts(num_nodes);
2412  std::vector<Real> magnitudes(num_nodes);
2413 
2414  for (int i=0; i<num_nodes; ++i)
2415  {
2416  real_parts[i] = local_soln[i].real();
2417  imag_parts[i] = local_soln[i].imag();
2418  magnitudes[i] = std::abs(local_soln[i]);
2419  }
2420 
2421  // Write the real, imaginary, and magnitude values to file.
2422  write_nodal_values(3 * variable_name_position + 1, real_parts, timestep);
2423  write_nodal_values(3 * variable_name_position + 2, imag_parts, timestep);
2424  write_nodal_values(3 * variable_name_position + 3, magnitudes, timestep);
2425 #endif
2426  }
2427 }
2428 
2429 
2430 
2432  const std::vector<std::pair<unsigned int, unsigned int>> & var_nums,
2433  int timestep,
2434  const std::vector<std::string> & output_names)
2435 {
2436  const MeshBase & mesh = es.get_mesh();
2437 
2438  // FIXME - half this code might be replaceable with a call to
2439  // EquationSystems::build_parallel_solution_vector()...
2440 
2441  for (auto [sys_num, var] : var_nums)
2442  {
2443  const System & sys = es.get_system(sys_num);
2444  const std::string & name = sys.variable_name(var);
2445 
2446  auto pos = std::find(output_names.begin(), output_names.end(), name);
2447 
2448  // Skip this name if it's not supposed to be output.
2449  if (pos == output_names.end())
2450  continue;
2451 
2452  // Compute the (zero-based) index which determines which
2453  // variable this will be as far as Nemesis is concerned. This
2454  // will be used below in the write_nodal_values() call.
2455  int variable_name_position =
2456  cast_int<int>(std::distance(output_names.begin(), pos));
2457 
2458  // Fill up a std::vector with the dofs for the current variable
2459  std::vector<numeric_index_type> required_indices(this->num_nodes);
2460 
2461  // Get the dof values required to write just our local part of
2462  // the solution vector.
2463  std::vector<Number> local_soln;
2464 
2465  const FEType type = sys.variable_type(var);
2466  if (type.family == SCALAR)
2467  {
2468  std::vector<numeric_index_type> scalar_indices;
2469  sys.get_dof_map().SCALAR_dof_indices(scalar_indices, var);
2470  for (int i=0; i<this->num_nodes; i++)
2471  required_indices[i] = scalar_indices[0];
2472  sys.current_local_solution->get(required_indices, local_soln);
2473  }
2474  else
2475  {
2476  // If we have DoFs at all nodes, e.g. for isoparametric
2477  // elements, this is easy:
2478  bool found_all_indices = true;
2479  for (int i=0; i<this->num_nodes; i++)
2480  {
2481  const Node & node = mesh.node_ref(this->exodus_node_num_to_libmesh[i]-1);
2482  if (node.n_comp(sys_num, var))
2483  required_indices[i] = node.dof_number(sys_num, var, 0);
2484  else
2485  {
2486  found_all_indices = false;
2487  break;
2488  }
2489  }
2490 
2491  if (found_all_indices)
2492  sys.current_local_solution->get(required_indices, local_soln);
2493  // Fine, we'll do it the hard way
2494  if (!found_all_indices)
2495  {
2496  local_soln.resize(num_nodes);
2497 
2498  const Variable & var_description = sys.variable(var);
2499  const DofMap & dof_map = sys.get_dof_map();
2500 
2502  std::vector<Number> elem_soln; // The finite element solution
2503  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
2504  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
2505 
2506  for (const auto & elem : mesh.active_local_element_ptr_range())
2507  if (var_description.active_on_subdomain(elem->subdomain_id()))
2508  {
2509  dof_map.dof_indices (elem, dof_indices, var);
2510  elem_soln.resize(dof_indices.size());
2511 
2512  for (auto i : index_range(dof_indices))
2513  elem_soln[i] = sys_soln(dof_indices[i]);
2514 
2515  FEInterface::nodal_soln (elem->dim(),
2516  type,
2517  elem,
2518  elem_soln,
2519  nodal_soln);
2520 
2521  // infinite elements should be skipped...
2522  if (!elem->infinite())
2523  for (auto n : elem->node_index_range())
2524  {
2525  const std::size_t exodus_num =
2526  libmesh_node_num_to_exodus[elem->node_id(n)];
2527  libmesh_assert_greater(exodus_num, 0);
2528  libmesh_assert_less(exodus_num-1, local_soln.size());
2529  local_soln[exodus_num-1] = nodal_soln[n];
2530  }
2531  }
2532  }
2533  }
2534 
2535 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
2536  // Call the ExodusII_IO_Helper function to write the data.
2537  write_nodal_values(variable_name_position + 1, local_soln, timestep);
2538 #else
2539  // We have the local (complex) values. Now extract the real,
2540  // imaginary, and magnitude values from them.
2541  std::vector<Real> real_parts(num_nodes);
2542  std::vector<Real> imag_parts(num_nodes);
2543  std::vector<Real> magnitudes(num_nodes);
2544 
2545  for (int i=0; i<num_nodes; ++i)
2546  {
2547  real_parts[i] = local_soln[i].real();
2548  imag_parts[i] = local_soln[i].imag();
2549  magnitudes[i] = std::abs(local_soln[i]);
2550  }
2551 
2552  // Write the real, imaginary, and magnitude values to file.
2553  write_nodal_values(3 * variable_name_position + 1, real_parts, timestep);
2554  write_nodal_values(3 * variable_name_position + 2, imag_parts, timestep);
2555  write_nodal_values(3 * variable_name_position + 3, magnitudes, timestep);
2556 #endif
2557  }
2558 }
2559 
2560 
2561 
2562 void
2564  const std::vector<std::set<subdomain_id_type>> & vars_active_subdomains)
2565 {
2566  // Quick return if there are no element variables to write
2567  if (names.size() == 0)
2568  return;
2569 
2570  // Quick return if we have already called this function
2572  return;
2573 
2574  // Be sure that variables in the file match what we are asking for
2575  if (num_elem_vars > 0)
2576  {
2577  this->check_existing_vars(ELEMENTAL, names, this->elem_var_names);
2578  return;
2579  }
2580 
2581  // Set the flag so we can skip this stuff on subsequent calls to
2582  // initialize_element_variables()
2583  _elem_vars_initialized = true;
2584 
2585  this->write_var_names(ELEMENTAL, names);
2586 
2587  // Create a truth table from global_elem_blk_ids and the information
2588  // in vars_active_subdomains. Create a truth table of
2589  // size global_elem_blk_ids.size() * names.size().
2590  std::vector<int> truth_tab(global_elem_blk_ids.size() * names.size());
2591  for (auto blk : index_range(global_elem_blk_ids))
2592  for (auto var : index_range(names))
2593  if (vars_active_subdomains[var].empty() ||
2594  vars_active_subdomains[var].count(cast_int<subdomain_id_type>(global_elem_blk_ids[blk])))
2595  truth_tab[names.size() * blk + var] = 1;
2596 
2597  // Write truth table to file.
2598  if (truth_tab.size())
2599  {
2600  ex_err = exII::ex_put_elem_var_tab(ex_id,
2601  cast_int<int>(global_elem_blk_ids.size()),
2602  cast_int<int>(names.size()),
2603  truth_tab.data());
2604  EX_CHECK_ERR(ex_err, "Error writing element truth table.");
2605  }
2606 }
2607 
2608 
2609 
2610 void
2612  const EquationSystems & es,
2613  const std::vector<std::pair<unsigned int, unsigned int>> &var_nums,
2614  int timestep,
2615  const std::vector<std::set<subdomain_id_type>> & vars_active_subdomains)
2616 {
2617  // For each variable in names,
2618  // For each subdomain in subdomain_map,
2619  // If this (subdomain, variable) combination is active
2620  // For each component in variable
2621  // Extract element values into local_soln (localize is a collective)
2622  // Write local_soln to file
2623  // Update var_ctr with number of vector components for variable
2624  //
2625  unsigned int var_ctr = 0;
2626  for (auto v : index_range(var_nums))
2627  {
2628  const unsigned int sys_num = var_nums[v].first;
2629  const unsigned int var = var_nums[v].second;
2630  const System & system = es.get_system(sys_num);
2631 
2632  // We need to check if the constant monomial is a scalar or a vector and set the number of
2633  // components as the mesh spatial dimension for the latter as per es.find_variable_numbers().
2634  // Even for the case where a variable is not active on any subdomain belonging to the
2635  // processor, we still need to know this number to update 'var_ctr'.
2636  const unsigned int n_comps =
2637  (system.variable_type(var) == FEType(CONSTANT, MONOMIAL_VEC)) ? mesh.spatial_dimension() : 1;
2638 
2639  // Get list of active subdomains for variable v
2640  const auto & active_subdomains = vars_active_subdomains[v];
2641 
2642  for (const int sbd_id_int : global_elem_blk_ids)
2643  {
2644  const subdomain_id_type sbd_id =
2645  cast_int<subdomain_id_type>(sbd_id_int);
2646  auto it = subdomain_map.find(sbd_id);
2647  const std::vector<dof_id_type> empty_vec;
2648  const std::vector<dof_id_type> & elem_ids =
2649  (it == subdomain_map.end()) ? empty_vec : it->second;
2650 
2651  // Possibly skip this (variable, subdomain) combination. Also, check that there is at
2652  // least one element on the subdomain... Indeed, it is possible to have zero elements,
2653  // e.g., when running "adaptivity_ex3" in parallel with the 'dimension=1' argument.
2654  if ((active_subdomains.empty() || active_subdomains.count(sbd_id)) && elem_ids.size())
2655  {
2656  std::vector<numeric_index_type> required_indices;
2657  required_indices.reserve(elem_ids.size());
2658 
2659  // The number of DOF components needs to be equal to the expected number so that we
2660  // know where to store data to correctly correspond to variable names - verify this by
2661  // accessing the n_comp method for the last element ID, which should return the same
2662  // value for all elements on a given subdomain, so we only need to check this once.
2663  libmesh_assert_equal_to(n_comps, mesh.elem_ref(elem_ids.back()).n_comp(sys_num, var));
2664 
2665  // Loop through the DOFs of the variable and write the values for it on each element.
2666  // The variable name should have been decomposed by es.find_variable_numbers().
2667  for (unsigned int comp = 0; comp < n_comps; ++comp)
2668  {
2669  for (const auto & id : elem_ids)
2670  required_indices.push_back(mesh.elem_ref(id).dof_number(sys_num, var, comp));
2671 
2672  std::vector<Number> local_soln;
2673  system.current_local_solution->get(required_indices, local_soln);
2674 
2675  // reset for the next component
2676  required_indices.clear();
2677 
2678  // It's possible that there's nothing for us to write:
2679  // we may not be responsible for any elements on the
2680  // current subdomain. We did still have to participate
2681  // in the localize() call above, however, since it is a
2682  // collective.
2683  if (local_soln.size())
2684  {
2685 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2686  int stride = write_complex_abs ? 3 : 2;
2687  std::vector<Real> local_soln_buffer(local_soln.size());
2688  std::transform(local_soln.begin(), local_soln.end(),
2689  local_soln_buffer.begin(), [](Number n) { return n.real(); });
2690  ex_err = exII::ex_put_elem_var(ex_id,
2691  timestep,
2692  static_cast<int>(stride*(var_ctr+comp)+1),
2693  static_cast<int>(sbd_id),
2694  static_cast<int>(local_soln.size()),
2695  MappedOutputVector(local_soln_buffer, _single_precision).data());
2696  EX_CHECK_ERR(ex_err, "Error writing element real values.");
2697 
2698  std::transform(local_soln.begin(), local_soln.end(),
2699  local_soln_buffer.begin(), [](Number n) { return n.imag(); });
2700  ex_err = exII::ex_put_elem_var(ex_id,
2701  timestep,
2702  static_cast<int>(stride*(var_ctr+comp)+2),
2703  static_cast<int>(sbd_id),
2704  static_cast<int>(local_soln.size()),
2705  MappedOutputVector(local_soln_buffer, _single_precision).data());
2706  EX_CHECK_ERR(ex_err, "Error writing element imaginary values.");
2707 
2708  if (write_complex_abs)
2709  {
2710  std::transform(local_soln.begin(), local_soln.end(),
2711  local_soln_buffer.begin(), [](Number n) { return std::abs(n); });
2712  ex_err = exII::ex_put_elem_var(ex_id,
2713  timestep,
2714  static_cast<int>(stride*(var_ctr+comp)+2),
2715  static_cast<int>(sbd_id),
2716  static_cast<int>(local_soln.size()),
2717  MappedOutputVector(local_soln_buffer, _single_precision).data());
2718  EX_CHECK_ERR(ex_err, "Error writing element magnitudes.");
2719  }
2720 #else // LIBMESH_USE_COMPLEX_NUMBERS
2721  ex_err = exII::ex_put_elem_var(ex_id,
2722  timestep,
2723  static_cast<int>(var_ctr+comp+1),
2724  static_cast<int>(sbd_id),
2725  static_cast<int>(local_soln.size()),
2726  MappedOutputVector(local_soln, _single_precision).data());
2727  EX_CHECK_ERR(ex_err, "Error writing element values.");
2728 #endif // LIBMESH_USE_COMPLEX_NUMBERS
2729  }
2730  } // end loop over vector components
2731  }
2732  } // end loop over active subdomains
2733 
2734  var_ctr += n_comps;
2735  } // end loop over vars
2736 
2737  this->update();
2738 }
2739 
2740 
2741 
2742 void Nemesis_IO_Helper::read_var_names_impl(const char * var_type,
2743  int & count,
2744  std::vector<std::string> & result)
2745 {
2746  // Most of what we need to do is the same as for Exodus
2747  this->ExodusII_IO_Helper::read_var_names_impl(var_type, count, result);
2748 
2749  // But with tests where we have more processors than elements,
2750  // Nemesis doesn't let us put variable names in files written by
2751  // processors owning nothing, but we may still *need* those
2752  // variable names on every processor, so let's sync them up...
2753 
2754  processor_id_type pid_broadcasting_names = this->processor_id();
2755  const std::size_t n_names = result.size();
2756  if (!n_names)
2757  pid_broadcasting_names = DofObject::invalid_processor_id;
2758 
2759  libmesh_assert(this->comm().semiverify
2760  (n_names ? nullptr : &n_names));
2761 
2762  this->comm().min(pid_broadcasting_names);
2763 
2764  if (pid_broadcasting_names != DofObject::invalid_processor_id)
2765  this->comm().broadcast(result, pid_broadcasting_names);
2766 }
2767 
2768 
2769 
2770 std::string Nemesis_IO_Helper::construct_nemesis_filename(std::string_view base_filename)
2771 {
2772  // Build a filename for this processor. This code is cut-n-pasted from the read function
2773  // and should probably be put into a separate function...
2774  std::ostringstream file_oss;
2775 
2776  // We have to be a little careful here: Nemesis left pads its file
2777  // numbers based on the number of processors, so for example on 10
2778  // processors, we'd have:
2779  // mesh.e.10.00
2780  // mesh.e.10.01
2781  // mesh.e.10.02
2782  // ...
2783  // mesh.e.10.09
2784 
2785  // And on 100 you'd have:
2786  // mesh.e.100.000
2787  // mesh.e.100.001
2788  // ...
2789  // mesh.e.128.099
2790 
2791  // Find the length of the highest processor ID
2792  file_oss << (this->n_processors());
2793  unsigned int field_width = cast_int<unsigned int>(file_oss.str().size());
2794 
2795  if (verbose)
2796  libMesh::out << "field_width=" << field_width << std::endl;
2797 
2798  file_oss.str(""); // reset the string stream
2799  file_oss << base_filename
2800  << '.' << this->n_processors()
2801  << '.' << std::setfill('0') << std::setw(field_width) << this->processor_id();
2802 
2803  // Return the resulting string
2804  return file_oss.str();
2805 }
2806 
2807 } // namespace libMesh
2808 
2809 #endif // #if defined(LIBMESH_HAVE_NEMESIS_API) && defined(LIBMESH_HAVE_EXODUS_API)
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
void write_var_names(ExodusVarType type, const std::vector< std::string > &names)
Wraps calls to exII::ex_put_var_names() and exII::ex_put_var_param().
FEFamily family
The type of finite element.
Definition: fe_type.h:221
This class facilitates inline conversion of an input data vector to a different precision level...
void allgather(const T &send_data, std::vector< T, A > &recv_data) const
This is the EquationSystems class.
const std::set< boundary_id_type > & get_side_boundary_ids() const
Nemesis_IO_Helper(const ParallelObject &parent, bool verbose=false, bool single_precision=false)
Constructor.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
std::map< subdomain_id_type, std::vector< dof_id_type > > subdomain_map
Map of subdomains to element numbers.
virtual ~Nemesis_IO_Helper()
Destructor.
std::vector< std::string > elem_var_names
void active_family_tree_by_side(std::vector< const Elem *> &family, unsigned int side, bool reset=true) const
Same as the active_family_tree() member, but only adds elements which are next to side...
Definition: elem.C:2167
A Node is like a Point, but with more information.
Definition: node.h:52
std::vector< int > num_global_node_counts
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
Definition: system.h:2458
void compute_element_maps()
This function computes element maps (really just packs vectors) which map the elements to internal an...
std::vector< int > global_nodeset_ids
Containers for reading global nodeset information.
virtual void initialize(std::string title, const MeshBase &mesh, bool use_discontinuous=false) override
Specialization of the initialize function from ExodusII_IO_Helper that also writes global initial dat...
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
void read_nodeset(int id)
Reading functions.
char ftype
The type of file to be written.
std::vector< int > global_sideset_ids
Containers for reading global sideset (boundary conditions) information.
void compute_num_global_nodesets(const MeshBase &pmesh)
This function uses global communication routines to determine the number of nodesets across the entir...
std::vector< int > node_cmap_ids
Vectors for storing the communication map parameters.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
void get_ss_param_global()
Fills: global_sideset_ids, num_global_side_counts, num_global_side_df_counts Call after: get_init_glo...
virtual void read_var_names_impl(const char *var_type, int &count, std::vector< std::string > &result)
read_var_names() dispatches to this function.
int num_external_nodes
The number of FEM nodes that reside on another processor but whose element partially resides on the c...
int nemesis_err_flag
Member data.
std::vector< int > exodus_elem_num_to_libmesh
std::vector< int > node_mape
Vector which stores external node IDs.
void sum(T &r) const
void compute_elem_communication_maps()
This function computes element communication maps (really just packs vectors) in preparation for writ...
void compute_num_global_elem_blocks(const MeshBase &pmesh)
This function uses global communication routines to determine the number of element blocks across the...
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
int num_node_cmaps
The number of nodal communication maps for this processor.
void put_elem_map(std::vector< int > &elem_mapi, std::vector< int > &elem_mapb)
Outputs IDs of internal and border elements.
const Parallel::Communicator & comm() const
std::vector< int > num_global_node_df_counts
std::vector< std::vector< int > > elem_cmap_side_ids
std::vector< std::vector< int > > node_cmap_proc_ids
void put_ns_param_global(std::vector< int > &global_nodeset_ids, std::vector< int > &num_global_node_counts, std::vector< int > &num_global_node_df_counts)
This function writes information about global node sets.
const ExodusII_IO_Helper::Conversion & get_conversion(const ElemType type) const
The libMesh namespace provides an interface to certain functionality in the library.
void put_eb_info_global(std::vector< int > &global_elem_blk_ids, std::vector< int > &global_elem_blk_cnts)
Writes global block information to the file .) global_elem_blk_ids - list of block IDs for all blocks...
std::vector< int > node_list
void compute_border_node_ids(const MeshBase &pmesh)
This function constructs the set of border node IDs present on the current mesh.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
bool write_complex_abs
By default, when complex numbers are enabled, for each variable we write out three values: the real p...
const T_sys & get_system(std::string_view name) const
void compute_communication_map_parameters()
This function determines the communication map parameters which will eventually be written to file...
void write_nodal_values(int var_id, const std::vector< Real > &values, int timestep)
Writes the vector of values to a nodal variable.
void message(std::string_view msg)
Prints the message defined in msg.
void check_existing_vars(ExodusVarType type, std::vector< std::string > &names, std::vector< std::string > &names_from_file)
When appending: during initialization, check that variable names in the file match those you attempt ...
Real distance(const Point &p)
virtual void write_elements(const MeshBase &mesh, bool use_discontinuous=false) override
This function is specialized to write the connectivity.
std::map< dof_id_type, dof_id_type > libmesh_node_num_to_exodus
int num_nodes_global
Global initial information.
void compute_node_communication_maps()
Compute the node communication maps (really just pack vectors) in preparation for writing them to fil...
uint8_t processor_id_type
Definition: id_types.h:104
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
std::vector< int > num_global_side_df_counts
This is the MeshBase class.
Definition: mesh_base.h:75
void close() noexcept
Closes the ExodusII mesh file.
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, sides, and ids for those sides.
void write_element_values(const MeshBase &mesh, const EquationSystems &es, const std::vector< std::pair< unsigned int, unsigned int >> &var_nums, int timestep, const std::vector< std::set< subdomain_id_type >> &vars_active_subdomains)
Writes the vector of elemental variable values, one variable and one subdomain at a time...
const std::set< boundary_id_type > & get_node_boundary_ids() const
void build_element_and_node_maps(const MeshBase &pmesh)
This function builds the libmesh -> exodus and exodus -> libmesh node and element maps...
void put_init_global(dof_id_type num_nodes_global, dof_id_type num_elems_global, unsigned num_elem_blks_global, unsigned num_node_sets_global, unsigned num_side_sets_global)
Writes global information including: .) global number of nodes .) global number of elems ...
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
processor_id_type n_processors() const
void build_node_list(std::vector< dof_id_type > &node_id_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of nodes and ids for those nodes.
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual void write_nodal_coordinates(const MeshBase &mesh, bool use_discontinuous=false) override
This function is specialized from ExodusII_IO_Helper to write only the nodal coordinates stored on th...
std::vector< int > num_global_side_counts
This class defines the notion of a variable in the system.
Definition: variable.h:50
void update()
Uses ex_update() to flush buffers to file.
void push_back_entry(const std::string &name)
Adds another name to the current data table.
void min(const T &r, T &o, Request &req) const
virtual unsigned int n_nodes() const =0
std::vector< int > global_elem_blk_ids
Read the global element block IDs and counts.
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
Definition: dof_object.h:493
virtual void initialize_element_variables(std::vector< std::string > names, const std::vector< std::set< subdomain_id_type >> &vars_active_subdomains) override
Override the Exodus Helper&#39;s implementation of this function so that it works correctly in parallel...
std::vector< std::vector< int > > elem_cmap_elem_ids
3 vectors of vectors for storing element communication IDs for this processor.
std::set< unsigned > border_elem_ids
A set of border elem IDs for this processor.
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
std::set< unsigned > internal_node_ids
A set of internal node IDs for this processor.
libmesh_assert(ctx)
std::vector< int > node_mapi
Vector which stores internal node IDs.
int num_proc
The number of processors for which the NEMESIS I file was created.
std::string construct_nemesis_filename(std::string_view base_filename)
Given base_filename, foo.e, constructs the Nemesis filename foo.e.X.Y, where X=n. ...
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2478
This is the ExodusII_IO_Helper class.
std::vector< std::vector< int > > elem_cmap_proc_ids
void put_cmap_params(std::vector< int > &node_cmap_ids, std::vector< int > &node_cmap_node_cnts, std::vector< int > &elem_cmap_ids, std::vector< int > &elem_cmap_elem_cnts)
Outputs initial information for communication maps.
std::set< unsigned > border_node_ids
The set which will eventually contain the IDs of "border nodes".
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:167
virtual dof_id_type parallel_n_nodes() const =0
std::set< int > nodes_attached_to_local_elems
libMesh numbered node ids attached to local elems.
void get_init_global()
Fills: num_nodes_global, num_elems_global, num_elem_blks_global, num_node_sets_global, num_side_sets_global Call after: read_and_store_header_info() Call before: Any other get_* function from this class.
An object whose state is distributed along a set of processors.
virtual dof_id_type parallel_n_elem() const =0
void compute_num_global_sidesets(const MeshBase &pmesh)
This function uses global communication routines to determine the number of sidesets across the entir...
std::map< dof_id_type, dof_id_type > libmesh_elem_num_to_exodus
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
int num_internal_elems
The number of internal FEM elements.
void compute_internal_and_border_elems_and_internal_nodes(const MeshBase &pmesh)
This function constructs the set of border and internal element IDs and internal node IDs present on ...
void compute_node_maps()
Compute the node maps (really just pack vectors) which map the nodes to internal, border...
void put_elem_cmap(std::vector< std::vector< int >> &elem_cmap_elem_ids, std::vector< std::vector< int >> &elem_cmap_side_ids, std::vector< std::vector< int >> &elem_cmap_proc_ids)
Writes information about elemental communication map.
std::vector< int > num_node_df_per_set
std::map< subdomain_id_type, unsigned > local_subdomain_counts
This map keeps track of the number of elements in each subdomain (block) for this processor...
std::map< int, std::vector< int > > block_id_to_elem_connectivity
This is the block connectivity, i.e.
char ** get_char_star_star()
Provide access to the underlying C data table.
virtual void write_nodesets(const MeshBase &mesh) override
Writes the nodesets for this processor.
std::map< unsigned, std::set< unsigned > >::iterator proc_nodes_touched_iterator
Typedef for an iterator into the data structure above.
std::unordered_set< const Node * > & node_set
Definition: mesh_tools.C:2276
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2508
void put_node_map(std::vector< int > &node_mapi, std::vector< int > &node_mapb, std::vector< int > &node_mape)
Outputs IDs of internal, border, and external nodes.
std::vector< int > elem_cmap_elem_cnts
int num_elem_cmaps
The number of elemental communication maps for this processor.
OStreamProxy out
std::vector< int > elem_cmap_ids
const MeshBase & get_mesh() const
static const bool value
Definition: xdr_io.C:54
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
void put_ss_param_global(std::vector< int > &global_sideset_ids, std::vector< int > &num_global_side_counts, std::vector< int > &num_global_side_df_counts)
This function writes information about global side sets.
virtual void write_sidesets(const MeshBase &mesh) override
Writes the sidesets for this processor.
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
int num_proc_in_file
The number of processors for which the NEMESIS I file stores information.
std::vector< int > node_cmap_node_cnts
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
std::map< unsigned, std::set< unsigned > > proc_nodes_touched_intersections
Another map to store sets of intersections with each other processor (other than ourself, of course).
std::vector< int > exodus_node_num_to_libmesh
void write_exodus_initialization_info(const MeshBase &pmesh, const std::string &title)
This function writes exodus-specific initialization information.
std::vector< int > num_nodes_per_set
void put_init_info(unsigned num_proc, unsigned num_proc_in_file, const char *ftype)
Writing functions.
void put_node_cmap(std::vector< std::vector< int >> &node_cmap_node_ids, std::vector< std::vector< int >> &node_cmap_proc_ids)
Outputs all of the nodal communication maps for this processor.
std::vector< int > elem_mapb
Vector which stores border element IDs.
std::vector< int > node_mapb
Vector which stores border node IDs.
virtual const Node * node_ptr(const dof_id_type i) const =0
std::vector< std::vector< int > > node_cmap_node_ids
2 vectors of vectors for storing the node communication IDs for this processor.
processor_id_type processor_id() const
std::map< unsigned, std::set< std::pair< unsigned, unsigned > > >::iterator proc_border_elem_sets_iterator
Typedef for an iterator into the data structure above.
This class is useful for managing anything that requires a char ** input/output in ExodusII file...
static void nodal_soln(const unsigned int dim, const FEType &fe_t, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, const bool add_p_level=true, const unsigned int vdim=1)
Build the nodal soln from the element soln.
Definition: fe_interface.C:625
const DofMap & get_dof_map() const
Definition: system.h:2374
processor_id_type processor_id() const
Definition: dof_object.h:905
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
std::vector< int > global_elem_blk_cnts
void put_loadbal_param(unsigned num_internal_nodes, unsigned num_border_nodes, unsigned num_external_nodes, unsigned num_internal_elems, unsigned num_border_elems, unsigned num_node_cmaps, unsigned num_elem_cmaps)
Writes load balance parameters, some of which are described below: .) num_internal_nodes - nodes "who...
std::vector< int > elem_mapi
Vector which stores internal element IDs.
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
std::map< unsigned, std::set< std::pair< unsigned, unsigned > > > proc_border_elem_sets
Map between processor ID and (element,side) pairs bordering that processor ID.
int num_border_nodes
The number of FEM nodes local to a processor but residing in an element which also has FEM nodes on o...
int num_internal_nodes
To be used with the Nemesis::ne_get_loadbal_param() routine.
std::set< unsigned > internal_elem_ids
A set of internal elem IDs for this processor.
uint8_t dof_id_type
Definition: id_types.h:67
virtual void read_var_names_impl(const char *var_type, int &count, std::vector< std::string > &result) override
read_var_names() dispatches to this function.
int num_border_elems
The number of border FEM elements.
void write_nodal_solution(const NumericVector< Number > &parallel_soln, const std::vector< std::string > &names, int timestep, const std::vector< std::string > &output_names)
Takes a parallel solution vector containing the node-major solution vector for all variables and outp...
void set_union(T &data, const unsigned int root_id) const
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.