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