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