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