Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // Local includes
20 : #include "libmesh/exodusII_io.h"
21 :
22 : #include "libmesh/boundary_info.h"
23 : #include "libmesh/cell_c0polyhedron.h"
24 : #include "libmesh/dof_map.h"
25 : #include "libmesh/dyna_io.h" // ElementDefinition for BEX
26 : #include "libmesh/enum_elem_type.h"
27 : #include "libmesh/elem.h"
28 : #include "libmesh/enum_to_string.h"
29 : #include "libmesh/equation_systems.h"
30 : #include "libmesh/exodusII_io_helper.h"
31 : #include "libmesh/face_c0polygon.h"
32 : #include "libmesh/fpe_disabler.h"
33 : #include "libmesh/int_range.h"
34 : #include "libmesh/libmesh_logging.h"
35 : #include "libmesh/mesh_base.h"
36 : #include "libmesh/mesh_communication.h"
37 : #include "libmesh/numeric_vector.h"
38 : #include "libmesh/parallel_mesh.h"
39 : #include "libmesh/parallel.h"
40 : #include "libmesh/system.h"
41 : #include "libmesh/utility.h"
42 :
43 : // TIMPI includes
44 : #include "timpi/parallel_sync.h"
45 :
46 : // C++ includes
47 : #include <cmath> // llround
48 : #include <cstring>
49 : #include <fstream>
50 : #include <map>
51 : #include <memory>
52 : #include <sstream>
53 :
54 : #ifdef LIBMESH_HAVE_EXODUS_API
55 : namespace
56 : {
57 : using namespace libMesh;
58 :
59 34608 : const std::vector<Real> & bex_constraint_vec(std::size_t i,
60 : const ExodusII_IO_Helper & helper)
61 : {
62 0 : std::size_t vec_offset = 0;
63 34608 : for (auto block_num : index_range(helper.bex_dense_constraint_vecs))
64 : {
65 0 : const auto & vecblock = helper.bex_dense_constraint_vecs[block_num];
66 0 : libmesh_assert_greater_equal(i, vec_offset);
67 34608 : if (i - vec_offset < vecblock.size())
68 34608 : return vecblock[i - vec_offset];
69 :
70 0 : vec_offset += vecblock.size();
71 : }
72 :
73 0 : libmesh_error_msg("Requested BEX coefficient vector " << i << " not found");
74 : }
75 :
76 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
77 : enum class SolutionOrdering
78 : {
79 : variable_major,
80 : geometry_major
81 : };
82 :
83 : std::size_t ordered_solution_index (const std::size_t var,
84 : const std::size_t geometry_index,
85 : const std::size_t num_vars,
86 : const std::size_t num_geometric_entries,
87 : const SolutionOrdering ordering)
88 : {
89 126 : if (ordering == SolutionOrdering::variable_major)
90 252 : return var * num_geometric_entries + geometry_index;
91 : else
92 0 : return geometry_index * num_vars + var;
93 : }
94 :
95 : std::vector<Real>
96 3 : complex_soln_components (const std::vector<Number> & soln,
97 : const std::size_t num_vars,
98 : const bool write_complex_abs,
99 : const SolutionOrdering source_ordering,
100 : const SolutionOrdering destination_ordering)
101 : {
102 : const std::size_t num_values = soln.size();
103 : libmesh_assert_greater(num_vars, 0);
104 : libmesh_assert_equal_to(num_values % num_vars, 0);
105 3 : const std::size_t num_geometric_entries = num_values / num_vars;
106 :
107 : // This will contain the real and imaginary parts and the magnitude
108 : // of the values in soln
109 3 : const std::size_t nco = write_complex_abs ? 3 : 2;
110 3 : const std::size_t num_complex_vars = nco * num_vars;
111 3 : std::vector<Real> complex_soln(nco * num_values);
112 :
113 10 : for (const auto var : make_range(num_vars))
114 70 : for (const auto geometry_index : make_range(num_geometric_entries))
115 : {
116 : const Number value =
117 : soln[ordered_solution_index
118 63 : (var, geometry_index, num_vars, num_geometric_entries, source_ordering)];
119 :
120 63 : const std::size_t complex_var = nco * var;
121 : complex_soln[ordered_solution_index
122 : (complex_var, geometry_index, num_complex_vars,
123 63 : num_geometric_entries, destination_ordering)] = value.real();
124 : complex_soln[ordered_solution_index
125 : (complex_var + 1, geometry_index, num_complex_vars,
126 63 : num_geometric_entries, destination_ordering)] = value.imag();
127 :
128 63 : if (write_complex_abs)
129 : complex_soln[ordered_solution_index
130 : (complex_var + 2, geometry_index, num_complex_vars,
131 126 : num_geometric_entries, destination_ordering)] = std::abs(value);
132 : }
133 :
134 3 : return complex_soln;
135 : }
136 : #endif // LIBMESH_USE_COMPLEX_NUMBERS
137 :
138 : }
139 : #endif
140 :
141 : namespace libMesh
142 : {
143 :
144 : // ------------------------------------------------------------
145 : // ExodusII_IO class members
146 31854 : ExodusII_IO::ExodusII_IO (MeshBase & mesh,
147 31854 : bool single_precision) :
148 : MeshInput<MeshBase> (mesh),
149 : MeshOutput<MeshBase> (mesh,
150 : /* is_parallel_format = */ false,
151 : /* serial_only_needed_on_proc_0 = */ true),
152 : ParallelObject(mesh),
153 : #ifdef LIBMESH_HAVE_EXODUS_API
154 0 : exio_helper(std::make_unique<ExodusII_IO_Helper>(*this, false, true, single_precision)),
155 29962 : _timestep(1),
156 29962 : _verbose(false),
157 29962 : _append(false),
158 : #endif
159 29962 : _allow_empty_variables(false),
160 29962 : _write_complex_abs(true),
161 29962 : _set_unique_ids_from_maps(false),
162 29962 : _disc_bex(false),
163 31854 : _allow_hi_order_elemental(false)
164 : {
165 : // if !LIBMESH_HAVE_EXODUS_API, we didn't use this
166 946 : libmesh_ignore(single_precision);
167 31854 : }
168 :
169 :
170 :
171 79 : ExodusII_IO::ExodusII_IO (const MeshBase & mesh,
172 79 : bool single_precision) :
173 : MeshInput<MeshBase> (),
174 : MeshOutput<MeshBase> (mesh,
175 : /* is_parallel_format = */ false,
176 : /* serial_only_needed_on_proc_0 = */ true),
177 : ParallelObject(mesh),
178 : #ifdef LIBMESH_HAVE_EXODUS_API
179 0 : exio_helper(std::make_unique<ExodusII_IO_Helper>(*this, false, true, single_precision)),
180 67 : _timestep(1),
181 67 : _verbose(false),
182 67 : _append(false),
183 : #endif
184 67 : _allow_empty_variables(false),
185 67 : _write_complex_abs(true),
186 67 : _set_unique_ids_from_maps(false),
187 67 : _disc_bex(false),
188 79 : _allow_hi_order_elemental(false)
189 : {
190 : // if !LIBMESH_HAVE_EXODUS_API, we didn't use this
191 6 : libmesh_ignore(single_precision);
192 79 : }
193 :
194 :
195 :
196 639 : int ExodusII_IO::get_exodus_version()
197 : {
198 : #ifdef LIBMESH_HAVE_EXODUS_API
199 639 : return ExodusII_IO_Helper::get_exodus_version();
200 : #else
201 : return 0;
202 : #endif
203 : }
204 :
205 :
206 :
207 142 : void ExodusII_IO::set_extra_integer_vars(const std::vector<std::string> & extra_integer_vars)
208 : {
209 142 : _extra_integer_vars = extra_integer_vars;
210 142 : }
211 :
212 70 : void ExodusII_IO::set_output_variables(const std::vector<std::string> & output_variables,
213 : bool allow_empty)
214 : {
215 70 : _output_variables = output_variables;
216 70 : _allow_empty_variables = allow_empty;
217 70 : }
218 :
219 :
220 :
221 284 : void ExodusII_IO::write_discontinuous_exodusII(const std::string & name,
222 : const EquationSystems & es,
223 : const std::set<std::string> * system_names)
224 : {
225 : // See write_element_data justification for serializing
226 : // Additionally, we need to serialize before building the subdomain name map
227 300 : MeshSerializer serialize(MeshInput<MeshBase>::mesh(), !MeshOutput<MeshBase>::_is_parallel_format, true);
228 :
229 24 : std::vector<std::string> solution_names;
230 16 : std::vector<Number> v;
231 :
232 284 : es.build_variable_names (solution_names, nullptr, system_names);
233 284 : es.build_discontinuous_solution_vector (v, system_names,
234 : nullptr, false, /* defaults */
235 284 : this->get_add_sides());
236 284 : this->write_nodal_data_discontinuous(name, v, solution_names);
237 284 : }
238 :
239 :
240 : #ifdef LIBMESH_HAVE_EXODUS_API
241 0 : void ExodusII_IO::write_timestep_discontinuous (const std::string &fname,
242 : const EquationSystems &es,
243 : const int timestep,
244 : const Real time,
245 : const std::set<std::string> * system_names)
246 : {
247 0 : _timestep = timestep;
248 0 : write_discontinuous_equation_systems (fname,es,system_names);
249 :
250 0 : if (MeshOutput<MeshBase>::mesh().processor_id())
251 0 : return;
252 :
253 0 : exio_helper->write_timestep(timestep, time);
254 : }
255 :
256 : #else
257 : void ExodusII_IO::write_timestep_discontinuous (const std::string & /* fname */,
258 : const EquationSystems & /* es */,
259 : const int /* timestep */,
260 : const Real /* time */,
261 : const std::set<std::string> * /*system_names*/)
262 : { libmesh_error(); }
263 : #endif
264 :
265 :
266 : // ------------------------------------------------------------
267 : // When the Exodus API is present...
268 : #ifdef LIBMESH_HAVE_EXODUS_API
269 :
270 32213 : ExodusII_IO::~ExodusII_IO ()
271 : {
272 31933 : exio_helper->close();
273 32213 : }
274 :
275 :
276 2482 : void ExodusII_IO::read (const std::string & fname)
277 : {
278 456 : LOG_SCOPE("read()", "ExodusII_IO");
279 :
280 : // Get a reference to the mesh we are reading
281 2482 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
282 :
283 : // Add extra integers into the mesh
284 304 : std::vector<unsigned int> extra_ids;
285 2624 : for (auto & name : _extra_integer_vars)
286 278 : extra_ids.push_back(mesh.add_elem_integer(name));
287 :
288 : // Clear any existing mesh data
289 2482 : mesh.clear();
290 :
291 : // Keep track of what kinds of elements this file contains
292 2482 : elems_of_dimension.clear();
293 2482 : elems_of_dimension.resize(4, false);
294 :
295 : // Open the exodus file in EX_READ mode
296 2482 : exio_helper->open(fname.c_str(), /*read_only=*/true);
297 :
298 : // Get header information from exodus file
299 2482 : exio_helper->read_and_store_header_info();
300 :
301 : // Read the QA records
302 2482 : exio_helper->read_qa_records();
303 :
304 : // Print header information
305 2482 : exio_helper->print_header();
306 :
307 : // Read nodes from the exodus file
308 2482 : exio_helper->read_nodes();
309 :
310 : // Reserve space for the nodes.
311 2634 : mesh.reserve_nodes(exio_helper->num_nodes);
312 :
313 : // Read the node number map from the Exodus file. This is
314 : // required if we want to preserve the numbering of nodes as it
315 : // exists in the Exodus file. If the Exodus file does not contain
316 : // a node_num_map, the identity map is returned by this call.
317 2482 : exio_helper->read_node_num_map();
318 :
319 : // Read the element number map from the Exodus file. This is
320 : // required if we want to preserve the numbering of elements as it
321 : // exists in the Exodus file. If the Exodus file does not contain
322 : // an elem_num_map, the identity map is returned by this call.
323 : //
324 : // We now do this before creating nodes, so if we have any spline
325 : // nodes that need a NodeElem attached we can give them unused elem
326 : // ids.
327 2482 : exio_helper->read_elem_num_map();
328 :
329 : // Read any Bezier Extraction coefficient vectors from the file,
330 : // such as might occur in an IsoGeometric Analysis (IGA) mesh.
331 2482 : exio_helper->read_bex_cv_blocks();
332 :
333 : // If we have Rational Bezier weights, we'll need to
334 : // store them.
335 152 : unsigned char weight_index = 0;
336 152 : const bool weights_exist = !exio_helper->w.empty();
337 :
338 : // If we have Bezier extraction coefficients, we'll need to put
339 : // NODEELEM elements on spline nodes, since our Rational Bezier
340 : // elements will be connected to nodes derived from those; nothing
341 : // else will be directly connected to the spline nodes.
342 152 : const bool bex_cv_exist = !exio_helper->bex_dense_constraint_vecs.empty();
343 :
344 : // If the user requests to set Node/Elem unique_ids based on the
345 : // node/elem_num_maps, then it is very likely that those unique_ids
346 : // will not be "unique" across the entire set of DofObjects (Nodes
347 : // and Elems), although they should of course still be unique within
348 : // the set of Nodes and the set of Elems on an individual basis. We
349 : // therefore need to relax some of the uniqueness checking that is
350 : // done in debug mode by setting the following flag on the Mesh.
351 2482 : if (_set_unique_ids_from_maps)
352 4 : mesh.allow_node_and_elem_unique_id_overlap(true);
353 :
354 : // Even if weights don't exist, we still use RATIONAL_BERNSTEIN for
355 : // Bezier-Bernstein BEX elements, we just use 1.0 as the weight on
356 : // every node.
357 2482 : if (bex_cv_exist)
358 : {
359 9 : const Real default_weight = 1.0;
360 : weight_index = cast_int<unsigned char>
361 85 : (mesh.add_node_datum<Real>("rational_weight", true,
362 : &default_weight));
363 0 : mesh.set_default_mapping_type(RATIONAL_BERNSTEIN_MAP);
364 0 : mesh.set_default_mapping_data(weight_index);
365 : }
366 :
367 304 : std::unordered_map<const Node *, Elem *> spline_nodeelem_ptrs;
368 :
369 : // Loop over the nodes, create Nodes with local processor_id 0.
370 328814 : for (auto i : make_range(exio_helper->num_nodes))
371 : {
372 : // Determine the libmesh node id implied by "i". The
373 : // get_libmesh_node_id() helper function expects a 1-based
374 : // Exodus node id, so we construct the "implied" Exodus node id
375 : // from "i" by adding 1.
376 326332 : auto libmesh_node_id = exio_helper->get_libmesh_node_id(/*exodus_node_id=*/i+1);
377 :
378 : // Catch the node that was added to the mesh
379 452992 : Node * added_node = mesh.add_point (Point(exio_helper->x[i], exio_helper->y[i], exio_helper->z[i]), libmesh_node_id);
380 :
381 : // Sanity check: throw an error if the Mesh assigned an ID to
382 : // the Node which does not match the libmesh_id we just determined.
383 326332 : libmesh_error_msg_if(added_node->id() != static_cast<unsigned>(libmesh_node_id),
384 : "Error! Mesh assigned node ID "
385 : << added_node->id()
386 : << " which is different from the (zero-based) Exodus ID "
387 : << libmesh_node_id
388 : << "!");
389 :
390 : // If the _set_unique_ids_from_maps flag is true, set the
391 : // unique_id for "node", otherwise do nothing.
392 326332 : exio_helper->conditionally_set_node_unique_id(mesh, added_node, i);
393 :
394 : // If we have a set of spline weights, these nodes are going to
395 : // be used as control points for Bezier elements, and we need
396 : // to attach a NodeElem to each to make sure it doesn't get
397 : // flagged as an unused node.
398 326332 : if (weights_exist)
399 : {
400 162 : const auto w = exio_helper->w[i];
401 162 : Point & p = *added_node;
402 0 : p /= w; // Exodus Bezier Extraction stores spline nodes in projective space
403 :
404 162 : added_node->set_extra_datum<Real>(weight_index, exio_helper->w[i]);
405 : }
406 :
407 326332 : if (bex_cv_exist)
408 : {
409 357 : std::unique_ptr<Elem> elem = Elem::build(NODEELEM);
410 :
411 : // Give the NodeElem ids at the end, so we can match any
412 : // existing ids in the file for other elements
413 : //
414 : // We don't set the unique_id for this NodeElem here even if
415 : // the user has set the _set_unique_ids_from_maps flag
416 : // because these NodeElems don't have entries in the
417 : // elem_num_map. Therefore, we just let the Mesh assign
418 : // whatever unique_id is "next" as the Elem is added to the
419 : // Mesh.
420 357 : elem->set_id() = exio_helper->end_elem_id() + i;
421 :
422 357 : elem->set_node(0, added_node);
423 357 : Elem * added_elem = mesh.add_elem(std::move(elem));
424 357 : spline_nodeelem_ptrs[added_node] = added_elem;
425 357 : }
426 : }
427 :
428 : // This assert is no longer valid if the nodes are not numbered
429 : // sequentially starting from 1 in the Exodus file.
430 : // libmesh_assert_equal_to (static_cast<unsigned int>(exio_helper->num_nodes), mesh.n_nodes());
431 :
432 : // Get information about all the element and edge blocks
433 2482 : exio_helper->read_block_info();
434 :
435 : // Reserve space for the elements. Account for any NodeElem that
436 : // have already been attached to spline control nodes.
437 2786 : mesh.reserve_elem(exio_helper->w.size() + exio_helper->num_elem);
438 :
439 : // Read variables for extra integer IDs
440 3005 : std::vector<std::map<dof_id_type, Real>> elem_ids(extra_ids.size());
441 : // We use the last time step to load the IDs
442 2482 : exio_helper->read_num_time_steps();
443 2482 : unsigned int last_step = exio_helper->num_time_steps;
444 2624 : for (auto i : index_range(extra_ids))
445 286 : exio_helper->read_elemental_var_values(_extra_integer_vars[i], last_step, elem_ids[i]);
446 :
447 : // Read in the element connectivity for each block.
448 152 : int nelem_last_block = 0;
449 :
450 : // If we're building Bezier elements from spline nodes, we need to
451 : // calculate those elements' local nodes on the fly, and we'll be
452 : // calculating them from constraint matrix columns, and we'll need
453 : // to make sure that the same node is found each time it's
454 : // calculated from multiple neighboring elements.
455 304 : std::map<std::vector<std::pair<dof_id_type, Real>>, Node *> local_nodes;
456 :
457 : // We'll set any spline NodeElem subdomain_id() values to exceed the
458 : // maximum of subdomain_id() values set via Exodus block ids.
459 2482 : subdomain_id_type max_subdomain_id = std::numeric_limits<subdomain_id_type>::min();
460 :
461 : // We've already added all the nodes explicitly specified in the
462 : // file, but if we have spline nodes we may need to add assembly
463 : // element nodes based on them. Add them contiguously so we're
464 : // compatible with any subsequent code paths (including our
465 : // ExodusII_IO::write()!) that don't support sparse ids.
466 2482 : dof_id_type n_nodes = mesh.n_nodes();
467 :
468 : // Loop over all the element blocks
469 9744 : for (int i=0; i<exio_helper->num_elem_blk; i++)
470 : {
471 : // Read the information for block i
472 7333 : exio_helper->read_elem_in_block (i);
473 : const subdomain_id_type subdomain_id =
474 7333 : restrict_int<subdomain_id_type>(exio_helper->get_block_id(i));
475 7333 : max_subdomain_id = std::max(max_subdomain_id, subdomain_id);
476 :
477 : // populate the map of names
478 7889 : std::string subdomain_name = exio_helper->get_block_name(i);
479 7333 : if (!subdomain_name.empty())
480 5524 : mesh.set_subdomain_name(subdomain_id, subdomain_name);
481 :
482 : // Set any relevant node/edge maps for this element
483 7402 : const std::string type_str (exio_helper->get_elem_type());
484 7402 : const auto & conv = exio_helper->get_conversion(type_str);
485 7333 : const bool is_c0polygon = (conv.libmesh_elem_type() == C0POLYGON);
486 7333 : const bool is_c0polyhedron = (conv.libmesh_elem_type() == C0POLYHEDRON);
487 556 : std::size_t c0polygon_connect_offset = 0;
488 556 : std::size_t c0polyhedron_connect_offset = 0;
489 :
490 : // Loop over all the faces in this block
491 7333 : int jmax = nelem_last_block+exio_helper->num_elem_this_blk;
492 311994 : for (int j=nelem_last_block; j<jmax; j++)
493 : {
494 304732 : const int elem_num = j - nelem_last_block;
495 279844 : std::unique_ptr<Elem> uelem;
496 :
497 304732 : if (is_c0polygon)
498 : {
499 12 : const int n_elem_nodes = exio_helper->elem_node_counts[elem_num];
500 12 : libmesh_error_msg_if(n_elem_nodes < 3,
501 : "Error: Exodus NSIDED block element "
502 : << elem_num
503 : << " has only " << n_elem_nodes << " nodes.");
504 :
505 13 : uelem = std::make_unique<C0Polygon>(cast_int<unsigned int>(n_elem_nodes));
506 : }
507 304720 : else if (is_c0polyhedron)
508 : {
509 36 : const int n_faces = exio_helper->elem_face_counts[elem_num];
510 36 : libmesh_error_msg_if(n_faces < 4,
511 : "Error: Exodus NFACED block element "
512 : << elem_num
513 : << " has only " << n_faces << " faces.");
514 :
515 : std::vector<std::shared_ptr<Polygon>>
516 42 : sides(cast_int<std::size_t>(n_faces));
517 :
518 324 : for (auto s : index_range(sides))
519 : {
520 24 : libmesh_assert_less(c0polyhedron_connect_offset,
521 : exio_helper->connect.size());
522 : const int exodus_face_id =
523 288 : exio_helper->connect[c0polyhedron_connect_offset++];
524 :
525 288 : libmesh_error_msg_if
526 : (exodus_face_id <= 0 ||
527 : exodus_face_id >
528 : cast_int<int>(exio_helper->c0polyhedron_face_connect.size()),
529 : "Error: Exodus NFACED block element "
530 : << elem_num
531 : << " references face ID " << exodus_face_id
532 : << ", but the file contains "
533 : << exio_helper->c0polyhedron_face_connect.size()
534 : << " faces.");
535 :
536 : const auto & face_nodes =
537 288 : exio_helper->c0polyhedron_face_connect[exodus_face_id - 1];
538 : auto side = std::make_shared<C0Polygon>
539 336 : (cast_int<unsigned int>(face_nodes.size()));
540 :
541 1584 : for (auto n : index_range(face_nodes))
542 : {
543 : const auto libmesh_node_id =
544 1404 : exio_helper->get_libmesh_node_id(face_nodes[n]);
545 1296 : side->set_node(n, mesh.node_ptr(libmesh_node_id));
546 : }
547 :
548 48 : sides[s] = std::move(side);
549 : }
550 :
551 36 : std::unique_ptr<Node> mid_elem_node;
552 39 : uelem = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
553 36 : if (mid_elem_node)
554 : {
555 39 : Node * added_node = mesh.add_node(std::move(mid_elem_node));
556 36 : if (added_node->id() >= n_nodes)
557 36 : n_nodes = added_node->id() + 1;
558 : }
559 30 : }
560 : else
561 584484 : uelem = Elem::build(conv.libmesh_elem_type());
562 :
563 : // Make sure that Exodus's number of nodes per Elem matches
564 : // the number of Nodes for this type of Elem. We only check
565 : // this for the first Elem in each block, since these values
566 : // are the same for every Elem in the block.
567 304732 : if (!is_c0polygon && !is_c0polyhedron && !elem_num)
568 7285 : libmesh_error_msg_if(exio_helper->num_nodes_per_elem != static_cast<int>(uelem->n_nodes()),
569 : "Error: Exodus file says "
570 : << exio_helper->num_nodes_per_elem
571 : << " nodes per Elem, but Elem type "
572 : << Utility::enum_to_string(uelem->type())
573 : << " has " << uelem->n_nodes() << " nodes.");
574 :
575 : // Assign the current subdomain to this Elem
576 304732 : uelem->subdomain_id() = subdomain_id;
577 :
578 : // Determine the libmesh elem id implied by "j". The
579 : // ExodusII_IO_Helper::get_libmesh_elem_id() helper function
580 : // expects a 1-based Exodus elem id, so we construct the
581 : // "implied" Exodus elem id from "j" by adding 1.
582 304732 : auto libmesh_elem_id = exio_helper->get_libmesh_elem_id(/*exodus_elem_id=*/j+1);
583 :
584 24888 : uelem->set_id(libmesh_elem_id);
585 :
586 : // Record that we have seen an element of dimension uelem->dim()
587 304732 : elems_of_dimension[uelem->dim()] = true;
588 :
589 : // Catch the Elem pointer that the Mesh throws back
590 354508 : Elem * elem = mesh.add_elem(std::move(uelem));
591 :
592 : // If the _set_unique_ids_from_maps flag is true, set the
593 : // unique_id for "elem", otherwise do nothing.
594 304732 : exio_helper->conditionally_set_elem_unique_id(mesh, elem, j);
595 :
596 : // If the Mesh assigned an ID different from the one we
597 : // tried to give it, we should probably error.
598 304732 : libmesh_error_msg_if(elem->id() != static_cast<unsigned>(libmesh_elem_id),
599 : "Error! Mesh assigned ID "
600 : << elem->id()
601 : << " which is different from the (zero-based) Exodus ID "
602 : << libmesh_elem_id
603 : << "!");
604 :
605 : // Assign extra integer IDs
606 305513 : for (auto & id : extra_ids)
607 : {
608 878 : const Real v = elem_ids[id][elem->id()];
609 :
610 852 : if (v == Real(-1))
611 : {
612 142 : elem->set_extra_integer(id, DofObject::invalid_id);
613 142 : continue;
614 : }
615 :
616 : // Ignore FE_INVALID here even if we've enabled FPEs; a
617 : // thrown exception is preferred over an FPE signal.
618 40 : FPEDisabler disable_fpes;
619 710 : const long long iv = std::llround(v);
620 :
621 : // Check if the real number is outside of the range we can
622 : // convert exactly
623 :
624 20 : long long max_representation = 1;
625 20 : max_representation = (max_representation << std::min(std::numeric_limits<Real>::digits,
626 20 : std::numeric_limits<double>::digits));
627 710 : libmesh_error_msg_if(iv > max_representation,
628 : "Error! An element integer value higher than "
629 : << max_representation
630 : << " was found! Exodus uses real numbers for storing element "
631 : " integers, which can only represent integers from 0 to "
632 : << max_representation
633 : << ".");
634 :
635 848 : libmesh_error_msg_if(iv < 0,
636 : "Error! An element integer value less than -1"
637 : << " was found! Exodus uses real numbers for storing element "
638 : " integers, which can only represent integers from 0 to "
639 : << max_representation
640 : << ".");
641 :
642 :
643 657 : elem->set_extra_integer(id, cast_int<dof_id_type>(iv));
644 : }
645 :
646 : // Set all the nodes for this element
647 : //
648 : // If we don't have any Bezier extraction operators, this
649 : // is easy: we've already built all our nodes and just need
650 : // to link to them.
651 304661 : if (is_c0polyhedron)
652 : {
653 : // Node pointers were already assigned while constructing
654 : // the polygonal sides above.
655 : }
656 304625 : else if (exio_helper->bex_cv_conn.empty())
657 : {
658 : const int n_nodes_this_elem =
659 329443 : is_c0polygon ?
660 12 : exio_helper->elem_node_counts[elem_num] :
661 24882 : exio_helper->num_nodes_per_elem;
662 :
663 1869560 : for (int k=0; k<n_nodes_this_elem; k++)
664 : {
665 : // Get index into this block's connectivity array
666 : std::size_t gi =
667 1565000 : is_c0polygon ?
668 60 : c0polygon_connect_offset++ :
669 1564940 : elem_num * exio_helper->num_nodes_per_elem + conv.get_node_map(k);
670 126623 : libmesh_assert_less(gi, exio_helper->connect.size());
671 :
672 : // Get the 1-based Exodus node id from the "connect" array
673 1565000 : auto exodus_node_id = exio_helper->connect[gi];
674 :
675 : // Convert this index to a libMesh Node id
676 1565000 : auto libmesh_node_id = exio_helper->get_libmesh_node_id(exodus_node_id);
677 :
678 : // Set the node pointer in the Elem
679 1565000 : elem->set_node(k, mesh.node_ptr(libmesh_node_id));
680 : }
681 : }
682 : else // We have Bezier Extraction data
683 : {
684 0 : auto & constraint_rows = mesh.set_constraint_rows();
685 :
686 : const DynaIO::ElementDefinition & dyna_elem_defn =
687 65 : DynaIO::find_elem_definition(elem->type(),
688 65 : elem->dim(),
689 65 : int(elem->default_order()));
690 :
691 : std::vector<std::vector<Real>>
692 67 : my_constraint_mat(exio_helper->bex_num_elem_cvs);
693 65 : for (auto spline_node_index :
694 1436 : make_range(exio_helper->bex_num_elem_cvs))
695 : {
696 1371 : my_constraint_mat[spline_node_index].resize(elem->n_nodes());
697 :
698 1371 : const auto & my_constraint_rows = exio_helper->bex_cv_conn[elem_num];
699 : const unsigned long elem_coef_vec_index =
700 1371 : my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based
701 1371 : const auto & my_vec = bex_constraint_vec(elem_coef_vec_index, *exio_helper);
702 0 : for (auto elem_node_index :
703 34608 : make_range(elem->n_nodes()))
704 : {
705 33237 : my_constraint_mat[spline_node_index][elem_node_index] =
706 33237 : my_vec[elem_node_index];
707 : }
708 :
709 : }
710 :
711 : // The tailing entries in each element's connectivity
712 : // vector are indices to Exodus constraint coefficient
713 : // rows.
714 :
715 : // Concatenating these rows gives a matrix with
716 : // all the constraints for the element nodes: each
717 : // column of that matrix is the constraint coefficients
718 : // for the node associated with that column (via the
719 : // Exodus numbering, not the libMesh numbering).
720 65 : const auto & my_constraint_rows = exio_helper->bex_cv_conn[elem_num];
721 :
722 0 : for (auto elem_node_index :
723 1424 : make_range(elem->n_nodes()))
724 : {
725 : // New finite element node data = dot product of
726 : // constraint matrix columns with spline node data.
727 : // Store each column's non-zero entries, along with
728 : // the global spline node indices, as a key to
729 : // identify shared finite element nodes.
730 0 : std::vector<std::pair<dof_id_type, Real>> key;
731 :
732 1359 : for (auto spline_node_index :
733 34596 : make_range(exio_helper->bex_num_elem_cvs))
734 : {
735 : // Pick out a row of the element constraint matrix
736 : const unsigned long elem_coef_vec_index =
737 33237 : my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based
738 :
739 : auto & coef_vec =
740 33237 : bex_constraint_vec(elem_coef_vec_index, *exio_helper);
741 :
742 : // Get coef from this node's column intersect that row
743 : const Real coef =
744 33237 : libmesh_vector_at(coef_vec, elem_node_index);
745 :
746 : // Get the libMesh node corresponding to that row
747 33237 : const int gi = elem_num * exio_helper->bex_num_elem_cvs + spline_node_index;
748 :
749 : // Get the 1-based Exodus node id from the "connect" array
750 33237 : auto exodus_node_id = exio_helper->connect[gi];
751 :
752 : // Convert this index to a libMesh Node id
753 33237 : auto libmesh_node_id = exio_helper->get_libmesh_node_id(exodus_node_id);
754 :
755 33237 : if (coef != 0) // Ignore irrelevant spline nodes
756 3661 : key.emplace_back(libmesh_node_id, coef);
757 : }
758 :
759 : // Have we already created this node? Connect it.
760 1359 : if (const auto local_node_it = local_nodes.find(key);
761 0 : local_node_it != local_nodes.end())
762 493 : elem->set_node(dyna_elem_defn.nodes[elem_node_index], local_node_it->second);
763 : // Have we not yet created this node? Construct it,
764 : // along with its weight and libMesh constraint row,
765 : // then connect it.
766 : else
767 : {
768 0 : Point p(0);
769 0 : Real w = 0;
770 0 : std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>> constraint_row;
771 :
772 2723 : for (auto [libmesh_spline_node_id, coef] : key)
773 : {
774 1857 : const Node & spline_node = mesh.node_ref(libmesh_spline_node_id);
775 :
776 1857 : p.add_scaled(spline_node, coef);
777 1857 : const Real spline_w = weights_exist ?
778 825 : spline_node.get_extra_datum<Real>(weight_index) : 1;
779 1857 : w += coef * spline_w;
780 :
781 : const Elem * nodeelem =
782 1857 : libmesh_map_find(spline_nodeelem_ptrs, &spline_node);
783 1857 : constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
784 : }
785 :
786 866 : Node *n = mesh.add_point(p, n_nodes++);
787 866 : if (weights_exist)
788 402 : n->set_extra_datum<Real>(weight_index, w);
789 :
790 : // If we're building disconnected Bezier
791 : // extraction elements then we don't want to
792 : // find the new nodes to reuse later; each
793 : // finite element node will connect to only one
794 : // element.
795 866 : if (!_disc_bex)
796 551 : local_nodes[key] = n;
797 866 : elem->set_node(dyna_elem_defn.nodes[elem_node_index], n);
798 :
799 866 : constraint_rows[n] = constraint_row;
800 : }
801 : }
802 65 : }
803 254956 : }
804 :
805 554 : libmesh_assert(!is_c0polygon ||
806 : c0polygon_connect_offset == exio_helper->connect.size());
807 554 : libmesh_assert(!is_c0polyhedron ||
808 : c0polyhedron_connect_offset == exio_helper->connect.size());
809 :
810 : // running sum of # of elements per block,
811 : // (should equal total number of elements in the end)
812 7262 : nelem_last_block += exio_helper->num_elem_this_blk;
813 : }
814 :
815 : // Now we know enough to fix any spline NodeElem subdomains
816 2411 : max_subdomain_id++;
817 2768 : for (auto p : spline_nodeelem_ptrs)
818 357 : p.second->subdomain_id() = max_subdomain_id;
819 :
820 : // Read in edge blocks, storing information in the BoundaryInfo object.
821 : // Edge blocks are treated as BCs.
822 2411 : exio_helper->read_edge_blocks(mesh);
823 :
824 : // Set the mesh dimension to the largest encountered for an element
825 12055 : for (unsigned char i=0; i!=4; ++i)
826 10244 : if (elems_of_dimension[i])
827 2843 : mesh.set_mesh_dimension(i);
828 :
829 : // Read in sideset information -- this is useful for applying boundary conditions
830 : {
831 : // Get basic information about all sidesets
832 2411 : exio_helper->read_sideset_info();
833 150 : int offset=0;
834 12397 : for (int i=0; i<exio_helper->num_side_sets; i++)
835 : {
836 : // Compute new offset
837 9986 : offset += (i > 0 ? exio_helper->num_sides_per_set[i-1] : 0);
838 9986 : exio_helper->read_sideset (i, offset);
839 :
840 10592 : std::string sideset_name = exio_helper->get_side_set_name(i);
841 9986 : if (!sideset_name.empty())
842 526 : mesh.get_boundary_info().sideset_name
843 8178 : (cast_int<boundary_id_type>(exio_helper->get_side_set_id(i)))
844 526 : = sideset_name;
845 : }
846 :
847 105191 : for (auto e : index_range(exio_helper->elem_list))
848 : {
849 : // Call helper function to get the libmesh Elem id for the
850 : // e'th entry in the current elem_list.
851 : dof_id_type libmesh_elem_id =
852 109980 : exio_helper->get_libmesh_elem_id(exio_helper->elem_list[e]);
853 :
854 : // Set any relevant node/edge maps for this element
855 102780 : Elem & elem = mesh.elem_ref(libmesh_elem_id);
856 :
857 102780 : const auto & conv = exio_helper->get_conversion(elem.type());
858 :
859 : // Map the zero-based Exodus side numbering to the libmesh side numbering
860 102780 : unsigned int raw_side_index = exio_helper->side_list[e]-1;
861 102780 : std::size_t side_index_offset = conv.get_shellface_index_offset();
862 :
863 102780 : if (raw_side_index < side_index_offset)
864 : {
865 : // We assume this is a "shell face"
866 3328 : int mapped_shellface = raw_side_index;
867 :
868 : // Check for errors
869 39936 : libmesh_error_msg_if(mapped_shellface < 0 || mapped_shellface >= 2,
870 : "Bad 0-based shellface id: "
871 : << mapped_shellface
872 : << " detected in Exodus file "
873 : << exio_helper->current_filename);
874 :
875 : // Add this (elem,shellface,id) triplet to the BoundaryInfo object.
876 39936 : mesh.get_boundary_info().add_shellface (libmesh_elem_id,
877 3328 : cast_int<unsigned short>(mapped_shellface),
878 43264 : cast_int<boundary_id_type>(exio_helper->id_list[e]));
879 : }
880 : else
881 : {
882 62844 : unsigned int side_index = static_cast<unsigned int>(raw_side_index - side_index_offset);
883 62844 : int mapped_side = conv.get_side_map(side_index);
884 :
885 : // Check for errors
886 62844 : libmesh_error_msg_if(mapped_side == ExodusII_IO_Helper::Conversion::invalid_id,
887 : "Invalid 1-based side id: "
888 : << side_index
889 : << " detected for "
890 : << Utility::enum_to_string(elem.type())
891 : << " in Exodus file "
892 : << exio_helper->current_filename);
893 :
894 62844 : libmesh_error_msg_if(mapped_side < 0 ||
895 : cast_int<unsigned int>(mapped_side) >= elem.n_sides(),
896 : "Bad 0-based side id: "
897 : << mapped_side
898 : << " detected for "
899 : << Utility::enum_to_string(elem.type())
900 : << " in Exodus file "
901 : << exio_helper->current_filename);
902 :
903 : // Add this (elem,side,id) triplet to the BoundaryInfo object.
904 62844 : mesh.get_boundary_info().add_side (libmesh_elem_id,
905 3872 : cast_int<unsigned short>(mapped_side),
906 66716 : cast_int<boundary_id_type>(exio_helper->id_list[e]));
907 : }
908 : } // end for (elem_list)
909 : } // end read sideset info
910 :
911 : // Read in elemset information and apply to Mesh elements if present
912 : {
913 2411 : exio_helper->read_elemset_info();
914 :
915 : // Mimic behavior of sideset case where we store all the set
916 : // information in a single array with offsets.
917 150 : int offset=0;
918 2553 : for (int i=0; i<exio_helper->num_elem_sets; i++)
919 : {
920 : // Compute new offset
921 142 : offset += (i > 0 ? exio_helper->num_elems_per_set[i-1] : 0);
922 142 : exio_helper->read_elemset (i, offset);
923 :
924 : // TODO: add support for elemset names
925 : // std::string elemset_name = exio_helper->get_elem_set_name(i);
926 : // if (!elemset_name.empty())
927 : // mesh.get_boundary_info().elemset_name(cast_int<boundary_id_type>(exio_helper->get_elem_set_id(i))) = elemset_name;
928 : }
929 :
930 : // Debugging: print the concatenated list of elemset ids
931 : // libMesh::out << "Concatenated list of elemset Elem ids (Exodus numbering):" << std::endl;
932 : // for (const auto & id : exio_helper->elemset_list)
933 : // libMesh::out << id << " ";
934 : // libMesh::out << std::endl;
935 :
936 : // Next we need to assign the elemset ids to the mesh using the
937 : // Elem's "extra_integers" support, if we have any.
938 2411 : if (exio_helper->num_elem_all_elemsets)
939 : {
940 : // Build map from Elem -> {elemsets}. This is needed only
941 : // temporarily to determine a unique set of elemset codes.
942 4 : std::map<Elem *, MeshBase::elemset_type> elem_to_elemsets;
943 639 : for (auto e : index_range(exio_helper->elemset_list))
944 : {
945 : // Call helper function to get the libmesh Elem id for the
946 : // e'th entry in the current elemset_list.
947 : dof_id_type libmesh_elem_id =
948 584 : exio_helper->get_libmesh_elem_id(exio_helper->elemset_list[e]);
949 :
950 : // Get a pointer to this Elem
951 568 : Elem * elem = mesh.elem_ptr(libmesh_elem_id);
952 :
953 : // Debugging:
954 : // libMesh::out << "Elem " << elem->id() << " is in elemset " << exio_helper->elemset_id_list[e] << std::endl;
955 :
956 : // Store elemset id in the map
957 568 : elem_to_elemsets[elem].insert(exio_helper->elemset_id_list[e]);
958 : }
959 :
960 : // Create a set of unique elemsets
961 4 : std::set<MeshBase::elemset_type> unique_elemsets;
962 497 : for (const auto & pr : elem_to_elemsets)
963 426 : unique_elemsets.insert(pr.second);
964 :
965 : // Debugging: print the unique elemsets
966 : // libMesh::out << "The set of unique elemsets which exist on the Mesh:" << std::endl;
967 : // for (const auto & s : unique_elemsets)
968 : // {
969 : // for (const auto & elemset_id : s)
970 : // libMesh::out << elemset_id << " ";
971 : // libMesh::out << std::endl;
972 : // }
973 :
974 : // Enumerate the unique_elemsets and tell the mesh about them
975 2 : dof_id_type code = 0;
976 284 : for (const auto & s : unique_elemsets)
977 420 : mesh.add_elemset_code(code++, s);
978 :
979 : // Sanity check: make sure that MeshBase::n_elemsets() reports
980 : // the expected value after calling MeshBase::add_elemset_code()
981 : // one or more times.
982 2 : libmesh_assert_msg(exio_helper->num_elem_sets == cast_int<int>(mesh.n_elemsets()),
983 : "Error: mesh.n_elemsets() is " << mesh.n_elemsets()
984 : << ", but mesh should have " << exio_helper->num_elem_sets << " elemsets.");
985 :
986 : // Create storage for the extra integer on all Elems. Elems which
987 : // are not in any set will use the default value of DofObject::invalid_id
988 : unsigned int elemset_index =
989 140 : mesh.add_elem_integer("elemset_code",
990 : /*allocate_data=*/true);
991 :
992 : // Store the appropriate extra_integer value on all Elems that need it.
993 497 : for (const auto & [elem, s] : elem_to_elemsets)
994 426 : elem->set_extra_integer(elemset_index, mesh.get_elemset_code(s));
995 : }
996 : } // done reading elemset info
997 :
998 : // Read nodeset info
999 : {
1000 : // This fills in the following fields of the helper for later use:
1001 : // nodeset_ids
1002 : // num_nodes_per_set
1003 : // num_node_df_per_set
1004 : // node_sets_node_index
1005 : // node_sets_dist_index
1006 : // node_sets_node_list
1007 : // node_sets_dist_fact
1008 2411 : exio_helper->read_all_nodesets();
1009 :
1010 10961 : for (int nodeset=0; nodeset<exio_helper->num_node_sets; nodeset++)
1011 : {
1012 : boundary_id_type nodeset_id =
1013 9108 : cast_int<boundary_id_type>(exio_helper->nodeset_ids[nodeset]);
1014 :
1015 9108 : std::string nodeset_name = exio_helper->get_node_set_name(nodeset);
1016 8550 : if (!nodeset_name.empty())
1017 8310 : mesh.get_boundary_info().nodeset_name(nodeset_id) = nodeset_name;
1018 :
1019 : // Get starting index of node ids for current nodeset.
1020 9108 : unsigned int offset = exio_helper->node_sets_node_index[nodeset];
1021 :
1022 110262 : for (int i=0; i<exio_helper->num_nodes_per_set[nodeset]; ++i)
1023 : {
1024 101154 : int exodus_node_id = exio_helper->node_sets_node_list[i + offset];
1025 101154 : auto libmesh_node_id = exio_helper->get_libmesh_node_id(exodus_node_id);
1026 101154 : mesh.get_boundary_info().add_node(libmesh_node_id, nodeset_id);
1027 : }
1028 : }
1029 : }
1030 :
1031 : #if LIBMESH_DIM < 3
1032 : libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
1033 : "Cannot open dimension "
1034 : << mesh.mesh_dimension()
1035 : << " mesh file when configured without "
1036 : << mesh.mesh_dimension()
1037 : << "D support.");
1038 : #endif
1039 4589 : }
1040 :
1041 :
1042 :
1043 : ExodusHeaderInfo
1044 213 : ExodusII_IO::read_header (const std::string & fname)
1045 : {
1046 : // We will need the Communicator of the Mesh we were created with.
1047 213 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1048 :
1049 : // Eventual return value
1050 6 : ExodusHeaderInfo header_info;
1051 :
1052 : // File I/O is done on processor 0, then broadcast to other procs
1053 219 : if (mesh.processor_id() == 0)
1054 : {
1055 : // Open the exodus file in EX_READ mode
1056 36 : exio_helper->open(fname.c_str(), /*read_only=*/true);
1057 :
1058 : // Get header information from exodus file without updating the
1059 : // Helper object's internal data structures.
1060 69 : header_info = exio_helper->read_header();
1061 :
1062 : // Close the file, we are now done with it. The goal is to keep the
1063 : // exio_helper object unchanged while calling this function,
1064 : // although it can't quite be marked "const" because we do have to
1065 : // actually open/close the file. This way, it should be possible to
1066 : // use the same ExodusII_IO object to read the headers of multiple
1067 : // different mesh files.
1068 36 : exio_helper->close();
1069 : }
1070 :
1071 : // Broadcast header_info to other procs before returning
1072 213 : header_info.broadcast(mesh.comm());
1073 :
1074 : // Return the information we read back to the user.
1075 213 : return header_info;
1076 : }
1077 :
1078 :
1079 :
1080 0 : void ExodusII_IO::verbose (bool set_verbosity)
1081 : {
1082 0 : _verbose = set_verbosity;
1083 :
1084 : // Set the verbose flag in the helper object as well.
1085 0 : exio_helper->verbose = _verbose;
1086 0 : }
1087 :
1088 :
1089 :
1090 0 : void ExodusII_IO::write_complex_magnitude (bool val)
1091 : {
1092 0 : _write_complex_abs = val;
1093 0 : }
1094 :
1095 284 : void ExodusII_IO::set_unique_ids_from_maps (bool val)
1096 : {
1097 284 : _set_unique_ids_from_maps = val;
1098 :
1099 : // Set this flag on the helper object as well. The helper needs to know about this
1100 : // flag, since it sometimes needs to construct libmesh Node ids from nodal connectivity
1101 : // arrays (see e.g. ExodusII_IO_Helper::read_edge_blocks()).
1102 284 : exio_helper->set_unique_ids_from_maps = val;
1103 284 : }
1104 :
1105 0 : void ExodusII_IO::use_mesh_dimension_instead_of_spatial_dimension(bool val)
1106 : {
1107 0 : exio_helper->use_mesh_dimension_instead_of_spatial_dimension(val);
1108 0 : }
1109 :
1110 :
1111 :
1112 0 : void ExodusII_IO::write_as_dimension(unsigned dim)
1113 : {
1114 0 : exio_helper->write_as_dimension(dim);
1115 0 : }
1116 :
1117 :
1118 :
1119 0 : void ExodusII_IO::set_coordinate_offset(Point p)
1120 : {
1121 : libmesh_warning("This method may be deprecated in the future");
1122 0 : exio_helper->set_coordinate_offset(p);
1123 0 : }
1124 :
1125 :
1126 :
1127 355 : void ExodusII_IO::append(bool val)
1128 : {
1129 355 : _append = val;
1130 355 : }
1131 :
1132 :
1133 :
1134 2698 : void ExodusII_IO::write_added_sides (bool val)
1135 : {
1136 76 : exio_helper->set_add_sides(val);
1137 2698 : }
1138 :
1139 :
1140 :
1141 49622 : bool ExodusII_IO::get_add_sides ()
1142 : {
1143 49622 : return exio_helper->get_add_sides();
1144 : }
1145 :
1146 :
1147 :
1148 :
1149 0 : const std::vector<Real> & ExodusII_IO::get_time_steps()
1150 : {
1151 0 : libmesh_error_msg_if
1152 : (!exio_helper->opened_for_reading,
1153 : "ERROR, ExodusII file must be opened for reading before calling ExodusII_IO::get_time_steps()!");
1154 :
1155 0 : exio_helper->read_time_steps();
1156 0 : return exio_helper->time_steps;
1157 : }
1158 :
1159 :
1160 :
1161 0 : int ExodusII_IO::get_num_time_steps()
1162 : {
1163 0 : libmesh_error_msg_if(!exio_helper->opened_for_reading && !exio_helper->opened_for_writing,
1164 : "ERROR, ExodusII file must be opened for reading or writing before calling ExodusII_IO::get_num_time_steps()!");
1165 :
1166 0 : exio_helper->read_num_time_steps();
1167 0 : return exio_helper->num_time_steps;
1168 : }
1169 :
1170 :
1171 :
1172 2994 : void ExodusII_IO::copy_nodal_solution(System & system,
1173 : std::string system_var_name,
1174 : std::string exodus_var_name,
1175 : unsigned int timestep)
1176 : {
1177 252 : LOG_SCOPE("copy_nodal_solution()", "ExodusII_IO");
1178 :
1179 2994 : const unsigned int var_num = system.variable_number(system_var_name);
1180 :
1181 2994 : const MeshBase & mesh = MeshInput<MeshBase>::mesh();
1182 :
1183 2994 : libmesh_error_msg_if(mesh.allow_renumbering(),
1184 : "ERROR, nodal data cannot be loaded if the mesh may be renumbered!");
1185 :
1186 : // With Exodus files we only open them on processor 0, so that's the
1187 : // where we have to do the data read too.
1188 2994 : if (system.comm().rank() == 0)
1189 : {
1190 516 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
1191 : "ERROR, ExodusII file must be opened for reading before copying a nodal solution!");
1192 :
1193 990 : exio_helper->read_nodal_var_values(exodus_var_name, timestep);
1194 : }
1195 :
1196 2994 : auto & node_var_value_map = exio_helper->nodal_var_values;
1197 :
1198 2994 : const bool serial_on_zero = mesh.is_serial_on_zero();
1199 :
1200 : // If our mesh isn't serial, then non-root processors need to
1201 : // request the data for their parts of the mesh and insert it
1202 : // themselves.
1203 2994 : if (!serial_on_zero)
1204 : {
1205 16 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> node_ids_to_request;
1206 2474 : if (this->processor_id() != 0)
1207 : {
1208 8 : std::vector<dof_id_type> node_ids;
1209 41425 : for (auto & node : mesh.local_node_ptr_range())
1210 39281 : node_ids.push_back(node->id());
1211 2140 : if (!node_ids.empty())
1212 3290 : node_ids_to_request[0] = std::move(node_ids);
1213 : }
1214 :
1215 : auto value_gather_functor =
1216 1639 : [& node_var_value_map]
1217 : (processor_id_type,
1218 : const std::vector<dof_id_type> & ids,
1219 37170 : std::vector<Real> & values)
1220 : {
1221 8 : const std::size_t query_size = ids.size();
1222 1647 : values.resize(query_size);
1223 38796 : for (std::size_t i=0; i != query_size; ++i)
1224 : {
1225 37162 : if (const auto it = node_var_value_map.find(ids[i]);
1226 13 : it != node_var_value_map.end())
1227 : {
1228 37162 : values[i] = it->second;
1229 37136 : node_var_value_map.erase(it);
1230 : }
1231 : else
1232 0 : values[i] = std::numeric_limits<Real>::quiet_NaN();
1233 : }
1234 4105 : };
1235 :
1236 : auto value_action_functor =
1237 1639 : [& node_var_value_map]
1238 : (processor_id_type,
1239 : const std::vector<dof_id_type> & ids,
1240 37144 : const std::vector<Real> & values)
1241 : {
1242 8 : const std::size_t query_size = ids.size();
1243 38796 : for (std::size_t i=0; i != query_size; ++i)
1244 37162 : if (!libmesh_isnan(values[i]))
1245 37162 : node_var_value_map[ids[i]] = values[i];
1246 2474 : };
1247 :
1248 8 : Real * value_ex = nullptr;
1249 : Parallel::pull_parallel_vector_data
1250 2466 : (system.comm(), node_ids_to_request, value_gather_functor,
1251 : value_action_functor, value_ex);
1252 : }
1253 :
1254 : // Everybody inserts the data they've received. If we're
1255 : // serial_on_zero then proc 0 inserts everybody's data and other
1256 : // procs have empty map ranges.
1257 78786 : for (auto p : exio_helper->nodal_var_values)
1258 : {
1259 6307 : dof_id_type i = p.first;
1260 75792 : const Node * node = MeshInput<MeshBase>::mesh().query_node_ptr(i);
1261 :
1262 75792 : if (node &&
1263 201965 : (serial_on_zero || node->processor_id() == system.processor_id()) &&
1264 75792 : node->n_comp(system.number(), var_num) > 0)
1265 : {
1266 75792 : dof_id_type dof_index = node->dof_number(system.number(), var_num, 0);
1267 :
1268 : // If the dof_index is local to this processor, set the value
1269 75792 : system.solution->set (dof_index, p.second);
1270 : }
1271 : }
1272 :
1273 2994 : system.solution->close();
1274 2994 : system.update();
1275 2994 : }
1276 :
1277 :
1278 :
1279 1687 : void ExodusII_IO::copy_elemental_solution(System & system,
1280 : std::string system_var_name,
1281 : std::string exodus_var_name,
1282 : unsigned int timestep)
1283 : {
1284 144 : LOG_SCOPE("copy_elemental_solution()", "ExodusII_IO");
1285 :
1286 1687 : const unsigned int var_num = system.variable_number(system_var_name);
1287 : // Assert that variable is an elemental data variable.
1288 : //
1289 : // NOTE: Currently, this reader is capable of reading only individual components of vector types,
1290 : // and each must be written out to its own scalar elemental data variable.
1291 1687 : const auto & var_type = system.variable_type(var_num);
1292 1687 : libmesh_error_msg_if(!EquationSystems::is_elemental_data_fe_type(var_type),
1293 : "Error! Trying to copy elemental solution into a variable that is not elemental data.");
1294 :
1295 1687 : const MeshBase & mesh = MeshInput<MeshBase>::mesh();
1296 48 : const DofMap & dof_map = system.get_dof_map();
1297 :
1298 1687 : libmesh_error_msg_if(mesh.allow_renumbering(),
1299 : "ERROR, elemental data cannot be loaded if the mesh may be renumbered!");
1300 :
1301 : // Map from element ID to elemental variable value. We need to use
1302 : // a map here rather than a vector (e.g. elem_var_values) since the
1303 : // libmesh element numbering can contain "holes". This is the case
1304 : // if we are reading elemental var values from an adaptively refined
1305 : // mesh that has not been sequentially renumbered.
1306 96 : std::map<dof_id_type, Real> elem_var_value_map;
1307 :
1308 : // With Exodus files we only open them on processor 0, so that's the
1309 : // where we have to do the data read too.
1310 1687 : if (system.comm().rank() == 0)
1311 : {
1312 271 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
1313 : "ERROR, ExodusII file must be opened for reading before copying an elemental solution!");
1314 :
1315 518 : exio_helper->read_elemental_var_values(exodus_var_name, timestep, elem_var_value_map);
1316 : }
1317 :
1318 1687 : const bool serial_on_zero = mesh.is_serial_on_zero();
1319 :
1320 : // If our mesh isn't serial, then non-root processors need to
1321 : // request the data for their parts of the mesh and insert it
1322 : // themselves.
1323 1687 : if (!serial_on_zero)
1324 : {
1325 28 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> elem_ids_to_request;
1326 507 : if (this->processor_id() != 0)
1327 : {
1328 14 : std::vector<dof_id_type> elem_ids;
1329 1274 : for (auto & elem : mesh.active_local_element_ptr_range())
1330 854 : elem_ids.push_back(elem->id());
1331 :
1332 413 : if (!elem_ids.empty())
1333 427 : elem_ids_to_request[0] = std::move(elem_ids);
1334 : }
1335 :
1336 : auto value_gather_functor =
1337 203 : [& elem_var_value_map]
1338 : (processor_id_type,
1339 : const std::vector<dof_id_type> & ids,
1340 497 : std::vector<Real> & values)
1341 : {
1342 14 : const std::size_t query_size = ids.size();
1343 217 : values.resize(query_size);
1344 672 : for (std::size_t i=0; i != query_size; ++i)
1345 : {
1346 483 : if (const auto it = elem_var_value_map.find(ids[i]);
1347 28 : it != elem_var_value_map.end())
1348 : {
1349 483 : values[i] = it->second;
1350 427 : elem_var_value_map.erase(it);
1351 : }
1352 : else
1353 0 : values[i] = std::numeric_limits<Real>::quiet_NaN();
1354 : }
1355 696 : };
1356 :
1357 : auto value_action_functor =
1358 203 : [& elem_var_value_map]
1359 : (processor_id_type,
1360 : const std::vector<dof_id_type> & ids,
1361 441 : const std::vector<Real> & values)
1362 : {
1363 14 : const std::size_t query_size = ids.size();
1364 672 : for (std::size_t i=0; i != query_size; ++i)
1365 483 : if (!libmesh_isnan(values[i]))
1366 483 : elem_var_value_map[ids[i]] = values[i];
1367 507 : };
1368 :
1369 14 : Real * value_ex = nullptr;
1370 : Parallel::pull_parallel_vector_data
1371 493 : (system.comm(), elem_ids_to_request, value_gather_functor,
1372 : value_action_functor, value_ex);
1373 : }
1374 :
1375 : std::map<dof_id_type, Real>::iterator
1376 48 : it = elem_var_value_map.begin(),
1377 48 : end = elem_var_value_map.end();
1378 :
1379 : // Everybody inserts the data they've received. If we're
1380 : // serial_on_zero then proc 0 inserts everybody's data and other
1381 : // procs have empty map ranges.
1382 4027 : for (; it!=end; ++it)
1383 : {
1384 2340 : const Elem * elem = mesh.query_elem_ptr(it->first);
1385 :
1386 2340 : if (elem && elem->n_comp(system.number(), var_num) > 0)
1387 : {
1388 2340 : dof_id_type dof_index = elem->dof_number(system.number(), var_num, 0);
1389 2340 : if (serial_on_zero || dof_map.local_index(dof_index ))
1390 2340 : system.solution->set (dof_index, it->second);
1391 : }
1392 : }
1393 :
1394 1687 : system.solution->close();
1395 1687 : system.update();
1396 1687 : }
1397 :
1398 0 : void ExodusII_IO::copy_scalar_solution(System & system,
1399 : std::vector<std::string> system_var_names,
1400 : std::vector<std::string> exodus_var_names,
1401 : unsigned int timestep)
1402 : {
1403 0 : LOG_SCOPE("copy_scalar_solution()", "ExodusII_IO");
1404 :
1405 0 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
1406 : "ERROR, ExodusII file must be opened for reading before copying a scalar solution!");
1407 :
1408 0 : libmesh_error_msg_if(system_var_names.size() != exodus_var_names.size(),
1409 : "ERROR, the number of system_var_names must match exodus_var_names.");
1410 :
1411 0 : std::vector<Real> values_from_exodus;
1412 0 : read_global_variable(exodus_var_names, timestep, values_from_exodus);
1413 :
1414 : #ifdef LIBMESH_HAVE_MPI
1415 0 : if (this->n_processors() > 1)
1416 : {
1417 0 : const Parallel::MessageTag tag = this->comm().get_unique_tag(1);
1418 0 : if (this->processor_id() == this->n_processors()-1)
1419 0 : this->comm().receive(0, values_from_exodus, tag);
1420 0 : if (this->processor_id() == 0)
1421 0 : this->comm().send(this->n_processors()-1, values_from_exodus, tag);
1422 0 : }
1423 : #endif
1424 :
1425 0 : if (system.processor_id() == (system.n_processors()-1))
1426 : {
1427 0 : const DofMap & dof_map = system.get_dof_map();
1428 :
1429 0 : for (auto i : index_range(system_var_names))
1430 : {
1431 0 : const unsigned int var_num = system.variable_scalar_number(system_var_names[i], 0);
1432 :
1433 0 : std::vector<dof_id_type> SCALAR_dofs;
1434 0 : dof_map.SCALAR_dof_indices(SCALAR_dofs, var_num);
1435 :
1436 0 : system.solution->set (SCALAR_dofs[0], values_from_exodus[i]);
1437 : }
1438 : }
1439 :
1440 0 : system.solution->close();
1441 0 : system.update();
1442 0 : }
1443 :
1444 0 : void ExodusII_IO::read_elemental_variable(std::string elemental_var_name,
1445 : unsigned int timestep,
1446 : std::map<unsigned int, Real> & unique_id_to_value_map)
1447 : {
1448 0 : LOG_SCOPE("read_elemental_variable()", "ExodusII_IO");
1449 :
1450 : // Note that this function MUST be called before renumbering
1451 0 : std::map<dof_id_type, Real> elem_var_value_map;
1452 :
1453 0 : exio_helper->read_elemental_var_values(elemental_var_name, timestep, elem_var_value_map);
1454 0 : for (auto & pr : elem_var_value_map)
1455 : {
1456 0 : const Elem * elem = MeshInput<MeshBase>::mesh().query_elem_ptr(pr.first);
1457 0 : unique_id_to_value_map.emplace(elem->top_parent()->unique_id(), pr.second);
1458 : }
1459 0 : }
1460 :
1461 0 : void ExodusII_IO::read_global_variable(std::vector<std::string> global_var_names,
1462 : unsigned int timestep,
1463 : std::vector<Real> & global_values)
1464 : {
1465 0 : LOG_SCOPE("read_global_variable()", "ExodusII_IO");
1466 :
1467 0 : std::size_t size = global_var_names.size();
1468 0 : libmesh_error_msg_if(size == 0, "ERROR, empty list of global variables to read from the Exodus file.");
1469 :
1470 : // read the values for all global variables
1471 0 : std::vector<Real> values_from_exodus;
1472 0 : exio_helper->read_var_names(ExodusII_IO_Helper::GLOBAL);
1473 0 : exio_helper->read_global_values(values_from_exodus, timestep);
1474 0 : std::vector<std::string> global_var_names_exodus = exio_helper->global_var_names;
1475 :
1476 0 : if (values_from_exodus.size() == 0)
1477 0 : return; // This will happen in parallel on procs that are not 0
1478 :
1479 0 : global_values.clear();
1480 0 : for (std::size_t i = 0; i != size; ++i)
1481 : {
1482 : // for each global variable in global_var_names, look the corresponding one in global_var_names_from_exodus
1483 : // and fill global_values accordingly
1484 0 : auto it = find(global_var_names_exodus.begin(), global_var_names_exodus.end(), global_var_names[i]);
1485 0 : if (it != global_var_names_exodus.end())
1486 0 : global_values.push_back(values_from_exodus[it - global_var_names_exodus.begin()]);
1487 : else
1488 0 : libmesh_error_msg("ERROR, Global variable " << global_var_names[i] << \
1489 : " not found in Exodus file.");
1490 : }
1491 :
1492 0 : }
1493 :
1494 423 : void ExodusII_IO::write_element_data (const EquationSystems & es)
1495 : {
1496 18 : LOG_SCOPE("write_element_data()", "ExodusII_IO");
1497 :
1498 : // Be sure the file has been opened for writing!
1499 435 : libmesh_error_msg_if(MeshOutput<MeshBase>::mesh().processor_id() == 0 && !exio_helper->opened_for_writing,
1500 : "ERROR, ExodusII file must be initialized before outputting element variables.");
1501 :
1502 : // This function currently only works on serialized meshes. We rely
1503 : // on having a reference to a non-const MeshBase object from our
1504 : // MeshInput parent class to construct a MeshSerializer object,
1505 : // similar to what is done in ExodusII_IO::write(). Note that
1506 : // calling ExodusII_IO::write_timestep() followed by
1507 : // ExodusII_IO::write_element_data() when the underlying Mesh is a
1508 : // DistributedMesh will result in an unnecessary additional
1509 : // serialization/re-parallelization step.
1510 : // The "true" specifies that we only need the mesh serialized to processor 0
1511 429 : MeshSerializer serialize(MeshInput<MeshBase>::mesh(), !MeshOutput<MeshBase>::_is_parallel_format, true);
1512 :
1513 : // To be (possibly) filled with a filtered list of variable names to output.
1514 18 : std::vector<std::string> names;
1515 :
1516 : // If _output_variables is populated, build_elemental_solution_vector() will filter this list to
1517 : // the variables that can be written as elemental data.
1518 423 : if (_output_variables.size() > 0)
1519 68 : names.assign(_output_variables.begin(), _output_variables.end());
1520 :
1521 : // If we pass in a list of names to "build_elemental_solution_vector()"
1522 : // it'll filter the variables coming back.
1523 12 : std::vector<Number> soln;
1524 423 : es.build_elemental_solution_vector(soln, names, _allow_hi_order_elemental);
1525 :
1526 : // Also, store the list of subdomains on which each variable is active
1527 18 : std::vector<std::set<subdomain_id_type>> vars_active_subdomains;
1528 423 : es.get_vars_active_subdomains(names, vars_active_subdomains);
1529 :
1530 423 : if (soln.empty()) // If there is nothing to write just return
1531 12 : return;
1532 :
1533 : // The data must ultimately be written block by block. This means that this data
1534 : // must be sorted appropriately.
1535 75 : if (MeshOutput<MeshBase>::mesh().processor_id())
1536 0 : return;
1537 :
1538 6 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1539 :
1540 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1541 :
1542 : std::vector<std::string> complex_names =
1543 3 : exio_helper->get_complex_names(names, _write_complex_abs);
1544 :
1545 : std::vector<std::set<subdomain_id_type>>
1546 : complex_vars_active_subdomains =
1547 : exio_helper->get_complex_vars_active_subdomains(vars_active_subdomains,
1548 3 : _write_complex_abs);
1549 3 : exio_helper->initialize_element_variables(complex_names, complex_vars_active_subdomains);
1550 :
1551 : const std::vector<Real> complex_soln =
1552 3 : complex_soln_components(soln, names.size(), _write_complex_abs,
1553 : SolutionOrdering::variable_major,
1554 3 : SolutionOrdering::variable_major);
1555 :
1556 3 : exio_helper->write_element_values(mesh, complex_soln, _timestep, complex_vars_active_subdomains);
1557 :
1558 : #else
1559 66 : exio_helper->initialize_element_variables(names, vars_active_subdomains);
1560 66 : exio_helper->write_element_values(mesh, soln, _timestep, vars_active_subdomains);
1561 : #endif
1562 798 : }
1563 :
1564 :
1565 :
1566 : void
1567 70 : ExodusII_IO::write_element_data_from_discontinuous_nodal_data
1568 : (const EquationSystems & es,
1569 : const std::set<std::string> * system_names,
1570 : const std::string & var_suffix)
1571 : {
1572 6 : LOG_SCOPE("write_element_data_from_discontinuous_nodal_data()", "ExodusII_IO");
1573 :
1574 : // Be sure that some other function has already opened the file and prepared it
1575 : // for writing. This is the same behavior as the write_element_data() function
1576 : // which we are trying to mimic.
1577 72 : libmesh_error_msg_if(MeshOutput<MeshBase>::mesh().processor_id() == 0 && !exio_helper->opened_for_writing,
1578 : "ERROR, ExodusII file must be initialized before outputting element variables.");
1579 :
1580 : // This function currently only works on serialized meshes. The
1581 : // "true" flag specifies that we only need the mesh serialized to
1582 : // processor 0
1583 : MeshSerializer serialize(MeshInput<MeshBase>::mesh(),
1584 70 : !MeshOutput<MeshBase>::_is_parallel_format,
1585 74 : true);
1586 :
1587 : // Note: in general we want to respect the contents of
1588 : // _output_variables, only building a solution vector with values
1589 : // from the requested variables. First build a list of all variable
1590 : // names, then throw out ones that aren't in _output_variables, if
1591 : // any.
1592 6 : std::vector<std::string> var_names;
1593 70 : es.build_variable_names (var_names, /*fetype=*/nullptr, system_names);
1594 :
1595 : // Get a subset of all variable names that can be written directly
1596 : // as elemental data. We treat those slightly differently since they
1597 : // truly only have a single value per Elem.
1598 6 : std::vector<std::string> elemental_var_names;
1599 70 : if (!_output_variables.empty())
1600 0 : elemental_var_names.assign(_output_variables.begin(), _output_variables.end());
1601 70 : es.build_elemental_data_variable_names(elemental_var_names, system_names);
1602 :
1603 : // Remove all names from var_names that are not in _output_variables.
1604 : // Note: This approach avoids errors when the user provides invalid
1605 : // variable names in _output_variables, as the code will not try to
1606 : // write a variable that doesn't exist.
1607 70 : if (!_output_variables.empty())
1608 : {
1609 : var_names.erase
1610 : (std::remove_if
1611 0 : (var_names.begin(),
1612 : var_names.end(),
1613 0 : [this](const std::string & name)
1614 0 : {return !std::count(_output_variables.begin(),
1615 : _output_variables.end(),
1616 0 : name);}),
1617 0 : var_names.end());
1618 : }
1619 :
1620 : // Build a solution vector, limiting the results to the variables in
1621 : // var_names and the Systems in system_names, and only computing values
1622 : // at the vertices.
1623 4 : std::vector<Number> v;
1624 : es.build_discontinuous_solution_vector
1625 70 : (v, system_names, &var_names, /*vertices_only=*/true);
1626 :
1627 : // Get active subdomains for each variable in var_names.
1628 6 : std::vector<std::set<subdomain_id_type>> vars_active_subdomains;
1629 70 : es.get_vars_active_subdomains(var_names, vars_active_subdomains);
1630 :
1631 : // Determine names of variables to write based on the number of
1632 : // nodes/vertices the elements in different subdomains have.
1633 4 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1634 4 : std::map<subdomain_id_type, unsigned int> subdomain_id_to_vertices_per_elem;
1635 970 : for (const auto & elem : mesh.active_element_ptr_range())
1636 : {
1637 : // Try to insert key/value pair into the map. If this returns
1638 : // false, check the returned iterator's value to make sure it
1639 : // matches. It shouldn't actually be possible for this to fail
1640 : // (since if the Mesh was like this it would have already
1641 : // failed) but it doesn't hurt to be on the safe side.
1642 : auto pr2 = subdomain_id_to_vertices_per_elem.emplace
1643 432 : (elem->subdomain_id(), elem->n_vertices());
1644 432 : libmesh_error_msg_if(!pr2.second && pr2.first->second != elem->n_vertices(),
1645 : "Elem with different number of vertices found.");
1646 66 : }
1647 :
1648 : // Determine "derived" variable names. These names are created by
1649 : // starting with the base variable name and appending the user's
1650 : // variable_suffix (default: "_elem_node_") followed by a node id.
1651 : //
1652 : // Not every derived variable will be active on every subdomain,
1653 : // even if the original variable _is_ active. Subdomains can have
1654 : // different geometric element types (with differing numbers of
1655 : // nodes), so some of the derived variable names will be inactive on
1656 : // those subdomains.
1657 : //
1658 : // Since we would otherwise generate the same name once per
1659 : // subdomain, we keep the list of names unique as we are creating
1660 : // it. We can't use a std::set for this because we don't want the
1661 : // variables names to be in a different order from the order
1662 : // they were written in the call to: build_discontinuous_solution_vector()
1663 : //
1664 : // The list of derived variable names includes one for each vertex,
1665 : // for higher-order elements we currently only write out vertex
1666 : // values, but this could be changed in the future without too much
1667 : // trouble.
1668 6 : std::vector<std::string> derived_var_names;
1669 :
1670 : // Keep track of mapping from derived_name to (orig_name, node_id)
1671 : // pair. We will use this later to determine whether a given
1672 : // variable is active on a given subdomain.
1673 : std::map<std::string, std::pair<std::string, unsigned int>>
1674 4 : derived_name_to_orig_name_and_node_id;
1675 :
1676 124 : for (const auto & pr : subdomain_id_to_vertices_per_elem)
1677 : {
1678 54 : const subdomain_id_type sbd_id = pr.first;
1679 : const unsigned int vertices_per_elem =
1680 54 : subdomain_id_to_vertices_per_elem[sbd_id];
1681 :
1682 58 : std::ostringstream oss;
1683 486 : for (unsigned int n=0; n<vertices_per_elem; ++n)
1684 3024 : for (const auto & orig_var_name : var_names)
1685 : {
1686 2592 : oss.str("");
1687 2592 : oss.clear();
1688 96 : oss << orig_var_name << var_suffix << n;
1689 192 : std::string derived_name = oss.str();
1690 :
1691 : // Only add this var name if it's not already in the list.
1692 2592 : if (!std::count(derived_var_names.begin(), derived_var_names.end(), derived_name))
1693 : {
1694 2592 : derived_var_names.push_back(derived_name);
1695 : // Add entry for derived_name -> (orig_name, node_id) mapping.
1696 2592 : derived_name_to_orig_name_and_node_id[derived_name] =
1697 2688 : std::make_pair(orig_var_name, n);
1698 : }
1699 : }
1700 50 : }
1701 :
1702 : // For each derived variable name, determine whether it is active
1703 : // based on how many nodes/vertices the elements in a given subdomain have,
1704 : // and whether they were active on the subdomain to begin with.
1705 : std::vector<std::set<subdomain_id_type>>
1706 76 : derived_vars_active_subdomains(derived_var_names.size());
1707 :
1708 : // A new data structure for keeping track of a list of variable names
1709 : // that are in the discontinuous solution vector on each subdomain. Used
1710 : // for indexing. Note: if a variable was inactive at the System level,
1711 : // an entry for it will still be in the discontinuous solution vector,
1712 : // but it will just have a value of zero. On the other hand, when we
1713 : // create the derived variable names some of them are "inactive" on
1714 : // different subdomains in the sense that they don't exist at all, i.e.
1715 : // there is no zero padding for them. We need to be able to distinguish
1716 : // between these two types in order to do the indexing into this vector
1717 : // correctly.
1718 : std::map<subdomain_id_type, std::vector<std::string>>
1719 4 : subdomain_to_var_names;
1720 :
1721 2662 : for (auto derived_var_id : index_range(derived_var_names))
1722 : {
1723 192 : const auto & derived_name = derived_var_names[derived_var_id];
1724 96 : const auto & [orig_name, node_id] =
1725 2592 : libmesh_map_find (derived_name_to_orig_name_and_node_id,
1726 : derived_name);
1727 :
1728 : // For each subdomain, determine whether the current variable
1729 : // should be active on that subdomain.
1730 5184 : for (const auto & pr : subdomain_id_to_vertices_per_elem)
1731 : {
1732 : // Convenience variables for the current subdomain and the
1733 : // number of nodes elements in this subdomain have.
1734 2592 : subdomain_id_type sbd_id = pr.first;
1735 : unsigned int vertices_per_elem_this_sbd =
1736 2592 : subdomain_id_to_vertices_per_elem[sbd_id];
1737 :
1738 : // Check whether variable orig_name was active on this
1739 : // subdomain to begin with by looking in the
1740 : // vars_active_subdomains container. We assume that the
1741 : // location of orig_name in the var_names vector matches its
1742 : // index in the vars_active_subdomains container.
1743 2592 : auto var_loc = std::find(var_names.begin(), var_names.end(), orig_name);
1744 2592 : libmesh_error_msg_if(var_loc == var_names.end(),
1745 : "Variable " << orig_name << " somehow not found in var_names array.");
1746 96 : auto var_id = std::distance(var_names.begin(), var_loc);
1747 :
1748 : // The derived_var will only be active if this subdomain has
1749 : // enough vertices for that to be the case.
1750 2592 : if (node_id < vertices_per_elem_this_sbd)
1751 : {
1752 : // Regardless of whether the original variable was not active on this subdomain,
1753 : // the discontinuous solution vector will have zero padding for it, and
1754 : // we will need to account for it. Therefore it should still be added to
1755 : // the subdomain_to_var_names data structure!
1756 2592 : subdomain_to_var_names[sbd_id].push_back(derived_name);
1757 :
1758 : // If the original variable was not active on the
1759 : // current subdomain, it should not be added to the
1760 : // derived_vars_active_subdomains data structure, since
1761 : // it will not be written to the Exodus file.
1762 :
1763 : // Determine if the original variable was active on the
1764 : // current subdomain.
1765 : bool orig_var_active =
1766 2688 : (vars_active_subdomains[var_id].empty() ||
1767 0 : vars_active_subdomains[var_id].count(sbd_id));
1768 :
1769 : // And only if it was, add it to the
1770 : // derived_vars_active_subdomains data structure.
1771 96 : if (orig_var_active)
1772 2592 : derived_vars_active_subdomains[derived_var_id].insert(sbd_id);
1773 : }
1774 : } // end loop over subdomain_id_to_vertices_per_elem
1775 : } // end loop over derived_var_names
1776 :
1777 : // At this point we've built the "true" list of derived names, but
1778 : // if there are any elemental data variables in this list, we now want
1779 : // to remove all but one copy of them from the derived_var_names list,
1780 : // and rename them in (but not remove them from) the subdomain_to_var_names
1781 : // list, and then update the derived_vars_active_subdomains containers
1782 : // before finally calling the Exodus helper functions.
1783 2662 : for (auto & derived_var_name : derived_var_names)
1784 : {
1785 : // Get the original name associated with this derived name.
1786 : const auto & name_and_id =
1787 2592 : libmesh_map_find (derived_name_to_orig_name_and_node_id,
1788 : derived_var_name);
1789 :
1790 : // Convenience variables for the map entry's contents.
1791 2592 : const std::string & orig_name = name_and_id.first;
1792 :
1793 : // Was the original name an elemental data variable?
1794 2592 : if (std::count(elemental_var_names.begin(),
1795 : elemental_var_names.end(),
1796 96 : orig_name))
1797 : {
1798 : // Rename this variable in the subdomain_to_var_names vectors.
1799 4320 : for (auto & pr : subdomain_to_var_names)
1800 : {
1801 : // Reference to ordered list of variable names on this subdomain.
1802 80 : auto & name_vec = pr.second;
1803 :
1804 : auto name_vec_it =
1805 2000 : std::find(name_vec.begin(),
1806 : name_vec.end(),
1807 160 : derived_var_name);
1808 :
1809 2160 : if (name_vec_it != name_vec.end())
1810 : {
1811 : // Actually rename it back to the orig_name, dropping
1812 : // the "_elem_corner_" stuff.
1813 80 : *name_vec_it = orig_name;
1814 : }
1815 : }
1816 :
1817 : // Finally, rename the variable in the derived_var_names vector itself.
1818 80 : derived_var_name = orig_name;
1819 : } // if (elemental data)
1820 : } // end loop over derived names
1821 :
1822 : // Now remove duplicate entries from derived_var_names after the first.
1823 : // Also update the derived_vars_active_subdomains container in a consistent way.
1824 : {
1825 6 : std::vector<std::string> derived_var_names_edited;
1826 6 : std::vector<std::set<subdomain_id_type>> derived_vars_active_subdomains_edited;
1827 74 : std::vector<unsigned int> found_first(elemental_var_names.size());
1828 :
1829 2662 : for (auto i : index_range(derived_var_names))
1830 : {
1831 192 : const auto & derived_var_name = derived_var_names[i];
1832 192 : const auto & active_set = derived_vars_active_subdomains[i];
1833 :
1834 : // Determine whether we will keep this derived variable name in
1835 : // the final container.
1836 96 : bool keep = true;
1837 15552 : for (auto j : index_range(elemental_var_names))
1838 12480 : if (derived_var_name == elemental_var_names[j])
1839 : {
1840 2240 : if (!found_first[j])
1841 270 : found_first[j] = 1;
1842 :
1843 : else
1844 70 : keep = false;
1845 : }
1846 :
1847 : // We also don't keep variables that are not active on any subdomains.
1848 : // Contrary to other uses of the var_active_subdomains container where
1849 : // the empty set means "all" subdomains, here it really means "none".
1850 2592 : if (active_set.empty())
1851 0 : keep = false;
1852 :
1853 2592 : if (keep)
1854 : {
1855 702 : derived_var_names_edited.push_back(derived_var_name);
1856 702 : derived_vars_active_subdomains_edited.push_back(active_set);
1857 : }
1858 : }
1859 :
1860 : // We built the filtered ranges, now swap them with the originals.
1861 2 : derived_var_names.swap(derived_var_names_edited);
1862 2 : derived_vars_active_subdomains.swap(derived_vars_active_subdomains_edited);
1863 66 : }
1864 :
1865 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1866 : // Build complex variable names "r_foo", "i_foo", "a_foo" and the lists of
1867 : // subdomains on which they are active.
1868 : auto complex_var_names =
1869 : exio_helper->get_complex_names(derived_var_names,
1870 0 : _write_complex_abs);
1871 : auto complex_vars_active_subdomains =
1872 : exio_helper->get_complex_vars_active_subdomains(derived_vars_active_subdomains,
1873 0 : _write_complex_abs);
1874 : auto complex_subdomain_to_var_names =
1875 : exio_helper->get_complex_subdomain_to_var_names(subdomain_to_var_names,
1876 0 : _write_complex_abs);
1877 :
1878 : // Make expanded version of vector "v" in which each entry in the
1879 : // original expands to an ("r_", "i_", "a_") triple.
1880 : // "nco" is the number of complex outputs, which depends on whether
1881 : // or not we are writing out the complex magnitudes.
1882 : std::vector<Real> complex_v;
1883 0 : int nco = _write_complex_abs ? 3 : 2;
1884 0 : complex_v.reserve(nco * v.size());
1885 0 : for (const auto & val : v)
1886 : {
1887 0 : complex_v.push_back(val.real());
1888 0 : complex_v.push_back(val.imag());
1889 0 : if (_write_complex_abs)
1890 0 : complex_v.push_back(std::abs(val));
1891 : }
1892 :
1893 : // Finally, initialize storage for the variables and write them to file.
1894 : exio_helper->initialize_element_variables
1895 0 : (complex_var_names, complex_vars_active_subdomains);
1896 : exio_helper->write_element_values_element_major
1897 0 : (mesh, complex_v, _timestep,
1898 : complex_vars_active_subdomains,
1899 : complex_var_names,
1900 : complex_subdomain_to_var_names);
1901 : #else
1902 :
1903 : // Call function which writes the derived variable names to the
1904 : // Exodus file.
1905 70 : exio_helper->initialize_element_variables(derived_var_names, derived_vars_active_subdomains);
1906 :
1907 : // ES::build_discontinuous_solution_vector() creates a vector with
1908 : // an element-major ordering, so call Helper::write_element_values()
1909 : // passing false for the last argument.
1910 : exio_helper->write_element_values_element_major
1911 70 : (mesh, v, _timestep,
1912 : derived_vars_active_subdomains,
1913 : derived_var_names,
1914 : subdomain_to_var_names);
1915 : #endif
1916 268 : }
1917 :
1918 :
1919 :
1920 48110 : void ExodusII_IO::write_nodal_data (const std::string & fname,
1921 : const std::vector<Number> & soln,
1922 : const std::vector<std::string> & names)
1923 : {
1924 2100 : LOG_SCOPE("write_nodal_data()", "ExodusII_IO");
1925 :
1926 2800 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1927 :
1928 2800 : int num_vars = cast_int<int>(names.size());
1929 48110 : dof_id_type num_nodes = mesh.n_nodes();
1930 :
1931 : // The names of the variables to be output
1932 2100 : std::vector<std::string> output_names;
1933 :
1934 48110 : if (_allow_empty_variables || !_output_variables.empty())
1935 0 : output_names = _output_variables;
1936 : else
1937 48110 : output_names = names;
1938 :
1939 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1940 : std::vector<std::string> complex_names =
1941 : exio_helper->get_complex_names(output_names,
1942 250 : _write_complex_abs);
1943 :
1944 : // Call helper function for opening/initializing data, giving it the
1945 : // complex variable names
1946 500 : this->write_nodal_data_common(fname, complex_names, /*continuous=*/true);
1947 : #else
1948 : // Call helper function for opening/initializing data
1949 92920 : this->write_nodal_data_common(fname, output_names, /*continuous=*/true);
1950 : #endif
1951 :
1952 49510 : if (mesh.processor_id())
1953 1400 : return;
1954 :
1955 : // This will count the number of variables actually output
1956 22234 : for (int c=0; c<num_vars; c++)
1957 : {
1958 15715 : std::stringstream name_to_find;
1959 :
1960 : std::vector<std::string>::iterator pos =
1961 15715 : std::find(output_names.begin(), output_names.end(), names[c]);
1962 14421 : if (pos == output_names.end())
1963 0 : continue;
1964 :
1965 : unsigned int variable_name_position =
1966 1294 : cast_int<unsigned int>(pos - output_names.begin());
1967 :
1968 : // Set up temporary vectors to be passed to Exodus to write the
1969 : // nodal values for a single variable at a time.
1970 : #ifdef LIBMESH_USE_REAL_NUMBERS
1971 2588 : std::vector<Number> cur_soln;
1972 :
1973 : // num_nodes is either exactly how much space we will need for
1974 : // each vector, or a safe upper bound for the amount of memory
1975 : // we will require when there are gaps in the numbering.
1976 13675 : cur_soln.reserve(num_nodes);
1977 : #else
1978 : std::vector<Real> real_parts;
1979 : std::vector<Real> imag_parts;
1980 : std::vector<Real> magnitudes;
1981 746 : real_parts.reserve(num_nodes);
1982 746 : imag_parts.reserve(num_nodes);
1983 746 : if (_write_complex_abs)
1984 746 : magnitudes.reserve(num_nodes);
1985 : #endif
1986 :
1987 : // There could be gaps in soln based on node numbering, but in
1988 : // serial the empty numbers are left empty.
1989 : // There could also be offsets in soln based on "fake" nodes
1990 : // inserted on each processor (because NumericVector indices
1991 : // have to be contiguous); the helper keeps track of those.
1992 : // We now copy the proper solution values contiguously into
1993 : // "cur_soln", removing the gaps.
1994 46342908 : for (const auto & node : mesh.node_ptr_range())
1995 : {
1996 : const dof_id_type idx =
1997 25252206 : (exio_helper->node_id_to_vec_id(node->id()))
1998 25252206 : * num_vars + c;
1999 : #ifdef LIBMESH_USE_REAL_NUMBERS
2000 25747343 : cur_soln.push_back(soln[idx]);
2001 : #else
2002 1599389 : real_parts.push_back(soln[idx].real());
2003 1599389 : imag_parts.push_back(soln[idx].imag());
2004 1599389 : if (_write_complex_abs)
2005 1599389 : magnitudes.push_back(std::abs(soln[idx]));
2006 : #endif
2007 11833 : }
2008 :
2009 : // If we're adding extra sides, we need to add their data too.
2010 : //
2011 : // Because soln was created from a parallel NumericVector, its
2012 : // numbering was contiguous on each processor; we need to use
2013 : // the same offsets here, and we need to loop through elements
2014 : // from earlier ranks first.
2015 14421 : if (exio_helper->get_add_sides())
2016 : {
2017 : std::vector<std::vector<const Elem *>>
2018 375 : elems_by_pid(mesh.n_processors());
2019 :
2020 3655 : for (const auto & elem : mesh.active_element_ptr_range())
2021 2070 : elems_by_pid[elem->processor_id()].push_back(elem);
2022 :
2023 2075 : for (auto p : index_range(elems_by_pid))
2024 : {
2025 : dof_id_type global_idx =
2026 1775 : exio_helper->added_node_offset_on(p) * num_vars + c;
2027 3455 : for (const Elem * elem : elems_by_pid[p])
2028 : {
2029 8880 : for (auto s : elem->side_index_range())
2030 : {
2031 7200 : if (EquationSystems::redundant_added_side(*elem,s))
2032 1944 : continue;
2033 :
2034 : const std::vector<unsigned int> side_nodes =
2035 5694 : elem->nodes_on_side(s);
2036 :
2037 38736 : for (auto n : index_range(side_nodes))
2038 : {
2039 2790 : libmesh_ignore(n);
2040 2790 : libmesh_assert_less(global_idx, soln.size());
2041 : #ifdef LIBMESH_USE_REAL_NUMBERS
2042 33480 : cur_soln.push_back(soln[global_idx]);
2043 : #else
2044 2790 : real_parts.push_back(soln[global_idx].real());
2045 2790 : imag_parts.push_back(soln[global_idx].imag());
2046 2790 : if (_write_complex_abs)
2047 2790 : magnitudes.push_back(std::abs(soln[global_idx]));
2048 : #endif
2049 33480 : global_idx += num_vars;
2050 : }
2051 : }
2052 : }
2053 : }
2054 250 : }
2055 :
2056 : // Finally, actually call the Exodus API to write to file.
2057 : #ifdef LIBMESH_USE_REAL_NUMBERS
2058 13675 : exio_helper->write_nodal_values(variable_name_position+1, cur_soln, _timestep);
2059 : #else
2060 746 : int nco = _write_complex_abs ? 3 : 2;
2061 746 : exio_helper->write_nodal_values(nco*variable_name_position+1, real_parts, _timestep);
2062 746 : exio_helper->write_nodal_values(nco*variable_name_position+2, imag_parts, _timestep);
2063 746 : if (_write_complex_abs)
2064 746 : exio_helper->write_nodal_values(3*variable_name_position+3, magnitudes, _timestep);
2065 : #endif
2066 :
2067 11833 : }
2068 45310 : }
2069 :
2070 :
2071 :
2072 :
2073 0 : void ExodusII_IO::write_information_records (const std::vector<std::string> & records)
2074 : {
2075 0 : if (MeshOutput<MeshBase>::mesh().processor_id())
2076 0 : return;
2077 :
2078 0 : libmesh_error_msg_if(!exio_helper->opened_for_writing,
2079 : "ERROR, ExodusII file must be initialized before outputting information records.");
2080 :
2081 0 : exio_helper->write_information_records(records);
2082 : }
2083 :
2084 :
2085 :
2086 0 : void ExodusII_IO::write_global_data (const std::vector<Number> & soln,
2087 : const std::vector<std::string> & names)
2088 : {
2089 0 : LOG_SCOPE("write_global_data()", "ExodusII_IO");
2090 :
2091 0 : if (MeshOutput<MeshBase>::mesh().processor_id())
2092 0 : return;
2093 :
2094 0 : libmesh_error_msg_if(!exio_helper->opened_for_writing,
2095 : "ERROR, ExodusII file must be initialized before outputting global variables.");
2096 :
2097 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2098 :
2099 : std::vector<std::string> complex_names =
2100 : exio_helper->get_complex_names(names,
2101 0 : _write_complex_abs);
2102 :
2103 0 : exio_helper->initialize_global_variables(complex_names);
2104 :
2105 : const std::vector<Real> complex_soln =
2106 0 : complex_soln_components(soln, names.size(), _write_complex_abs,
2107 : SolutionOrdering::variable_major,
2108 0 : SolutionOrdering::variable_major);
2109 :
2110 0 : exio_helper->write_global_values(complex_soln, _timestep);
2111 :
2112 : #else
2113 0 : exio_helper->initialize_global_variables(names);
2114 0 : exio_helper->write_global_values(soln, _timestep);
2115 : #endif
2116 0 : }
2117 :
2118 :
2119 :
2120 48110 : void ExodusII_IO::write_timestep (const std::string & fname,
2121 : const EquationSystems & es,
2122 : const int timestep,
2123 : const Real time,
2124 : const std::set<std::string> * system_names)
2125 : {
2126 48110 : _timestep = timestep;
2127 48110 : MeshOutput<MeshBase>::write_equation_systems(fname,es,system_names);
2128 :
2129 49510 : if (MeshOutput<MeshBase>::mesh().processor_id())
2130 700 : return;
2131 :
2132 7813 : exio_helper->write_timestep(timestep, time);
2133 : }
2134 :
2135 :
2136 13030 : void ExodusII_IO::write_equation_systems (const std::string & fname,
2137 : const EquationSystems & es,
2138 : const std::set<std::string> * system_names)
2139 : {
2140 13030 : write_timestep(fname, es, 1, 0, system_names);
2141 13030 : }
2142 :
2143 :
2144 0 : void ExodusII_IO::write_elemsets()
2145 : {
2146 0 : libmesh_error_msg_if(!exio_helper->opened_for_writing,
2147 : "ERROR, ExodusII file must be opened for writing "
2148 : "before calling ExodusII_IO::write_elemsets()!");
2149 :
2150 0 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
2151 0 : exio_helper->write_elemsets(mesh);
2152 0 : }
2153 :
2154 : void
2155 71 : ExodusII_IO::
2156 : write_sideset_data(int timestep,
2157 : const std::vector<std::string> & var_names,
2158 : const std::vector<std::set<boundary_id_type>> & side_ids,
2159 : const std::vector<std::map<BoundaryInfo::BCTuple, Real>> & bc_vals)
2160 : {
2161 71 : libmesh_error_msg_if(!exio_helper->opened_for_writing,
2162 : "ERROR, ExodusII file must be opened for writing "
2163 : "before calling ExodusII_IO::write_sideset_data()!");
2164 :
2165 4 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
2166 71 : exio_helper->write_sideset_data(mesh, timestep, var_names, side_ids, bc_vals);
2167 71 : }
2168 :
2169 :
2170 :
2171 : void
2172 71 : ExodusII_IO::
2173 : read_sideset_data(int timestep,
2174 : std::vector<std::string> & var_names,
2175 : std::vector<std::set<boundary_id_type>> & side_ids,
2176 : std::vector<std::map<BoundaryInfo::BCTuple, Real>> & bc_vals)
2177 : {
2178 71 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
2179 : "ERROR, ExodusII file must be opened for reading "
2180 : "before calling ExodusII_IO::read_sideset_data()!");
2181 :
2182 4 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
2183 71 : exio_helper->read_sideset_data(mesh, timestep, var_names, side_ids, bc_vals);
2184 71 : }
2185 :
2186 :
2187 :
2188 : void
2189 142 : ExodusII_IO::
2190 : get_sideset_data_indices (std::map<BoundaryInfo::BCTuple, unsigned int> & bc_array_indices)
2191 :
2192 : {
2193 142 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
2194 : "ERROR, ExodusII file must be opened for reading "
2195 : "before calling ExodusII_IO::get_sideset_data_indices()!");
2196 :
2197 8 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
2198 142 : exio_helper->get_sideset_data_indices(mesh, bc_array_indices);
2199 142 : }
2200 :
2201 : void
2202 142 : ExodusII_IO::
2203 : get_nodeset_data_indices (std::map<BoundaryInfo::NodeBCTuple, unsigned int> & bc_array_indices)
2204 : {
2205 142 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
2206 : "ERROR, ExodusII file must be opened for reading "
2207 : "before calling ExodusII_IO::get_nodeset_data_indices()!");
2208 :
2209 142 : exio_helper->get_nodeset_data_indices(bc_array_indices);
2210 142 : }
2211 :
2212 : void
2213 71 : ExodusII_IO::
2214 : write_nodeset_data (int timestep,
2215 : const std::vector<std::string> & var_names,
2216 : const std::vector<std::set<boundary_id_type>> & node_boundary_ids,
2217 : const std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> & bc_vals)
2218 : {
2219 71 : libmesh_error_msg_if(!exio_helper->opened_for_writing,
2220 : "ERROR, ExodusII file must be opened for writing "
2221 : "before calling ExodusII_IO::write_nodeset_data()!");
2222 :
2223 71 : exio_helper->write_nodeset_data(timestep, var_names, node_boundary_ids, bc_vals);
2224 71 : }
2225 :
2226 :
2227 :
2228 : void
2229 71 : ExodusII_IO::
2230 : read_nodeset_data (int timestep,
2231 : std::vector<std::string> & var_names,
2232 : std::vector<std::set<boundary_id_type>> & node_boundary_ids,
2233 : std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> & bc_vals)
2234 : {
2235 71 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
2236 : "ERROR, ExodusII file must be opened for reading "
2237 : "before calling ExodusII_IO::read_nodeset_data()!");
2238 :
2239 71 : exio_helper->read_nodeset_data(timestep, var_names, node_boundary_ids, bc_vals);
2240 71 : }
2241 :
2242 : void
2243 71 : ExodusII_IO::
2244 : write_elemset_data (int timestep,
2245 : const std::vector<std::string> & var_names,
2246 : const std::vector<std::set<elemset_id_type>> & elemset_ids_in,
2247 : const std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> & elemset_vals)
2248 : {
2249 71 : libmesh_error_msg_if(!exio_helper->opened_for_writing,
2250 : "ERROR, ExodusII file must be opened for writing "
2251 : "before calling ExodusII_IO::write_elemset_data()!");
2252 :
2253 71 : exio_helper->write_elemset_data(timestep, var_names, elemset_ids_in, elemset_vals);
2254 71 : }
2255 :
2256 :
2257 :
2258 : void
2259 71 : ExodusII_IO::
2260 : read_elemset_data (int timestep,
2261 : std::vector<std::string> & var_names,
2262 : std::vector<std::set<elemset_id_type>> & elemset_ids_in,
2263 : std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> & elemset_vals)
2264 : {
2265 71 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
2266 : "ERROR, ExodusII file must be opened for reading "
2267 : "before calling ExodusII_IO::read_elemset_data()!");
2268 :
2269 71 : exio_helper->read_elemset_data(timestep, var_names, elemset_ids_in, elemset_vals);
2270 71 : }
2271 :
2272 : void
2273 71 : ExodusII_IO::get_elemset_data_indices (std::map<std::pair<dof_id_type, elemset_id_type>, unsigned int> & elemset_array_indices)
2274 : {
2275 71 : libmesh_error_msg_if(!exio_helper->opened_for_reading,
2276 : "ERROR, ExodusII file must be opened for reading "
2277 : "before calling ExodusII_IO::get_elemset_data_indices()!");
2278 :
2279 71 : exio_helper->get_elemset_data_indices(elemset_array_indices);
2280 71 : }
2281 :
2282 :
2283 3301 : void ExodusII_IO::write (const std::string & fname)
2284 : {
2285 300 : LOG_SCOPE("write()", "ExodusII_IO");
2286 :
2287 200 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
2288 :
2289 : // We may need to gather a DistributedMesh to output it, making that
2290 : // const qualifier in our constructor a dirty lie
2291 : // The "true" specifies that we only need the mesh serialized to processor 0
2292 : MeshSerializer serialize
2293 : (const_cast<MeshBase &>(mesh),
2294 3501 : !MeshOutput<MeshBase>::_is_parallel_format, true);
2295 :
2296 100 : libmesh_assert( !exio_helper->opened_for_writing );
2297 :
2298 : // If the user has set the append flag here, it doesn't really make
2299 : // sense: the intent of this function is to write a Mesh with no
2300 : // data, while "appending" is really intended to add data to an
2301 : // existing file. If we're verbose, print a message to this effect.
2302 100 : if (_append && _verbose)
2303 : libmesh_warning("Warning: Appending in ExodusII_IO::write() does not make sense.\n"
2304 : "Creating a new file instead!");
2305 :
2306 6502 : exio_helper->create(fname);
2307 6502 : exio_helper->initialize(fname,mesh);
2308 3301 : exio_helper->write_nodal_coordinates(mesh);
2309 3301 : exio_helper->write_elements(mesh);
2310 3301 : exio_helper->write_sidesets(mesh);
2311 3301 : exio_helper->write_nodesets(mesh);
2312 3301 : exio_helper->write_elemsets(mesh);
2313 :
2314 3301 : if ((mesh.get_boundary_info().n_edge_conds() > 0) && _verbose)
2315 : libmesh_warning("Warning: Mesh contains edge boundary IDs, but these "
2316 : "are not supported by the ExodusII format.");
2317 3301 : }
2318 :
2319 :
2320 :
2321 1491 : void ExodusII_IO::write_nodal_data_discontinuous (const std::string & fname,
2322 : const std::vector<Number> & soln,
2323 : const std::vector<std::string> & names)
2324 : {
2325 63 : LOG_SCOPE("write_nodal_data_discontinuous()", "ExodusII_IO");
2326 :
2327 84 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
2328 :
2329 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2330 :
2331 : std::vector<std::string> complex_names =
2332 : exio_helper->get_complex_names(names,
2333 21 : _write_complex_abs);
2334 :
2335 : // Call helper function for opening/initializing data, giving it the
2336 : // complex variable names
2337 42 : this->write_nodal_data_common(fname, complex_names, /*continuous=*/false);
2338 : #else
2339 : // Call helper function for opening/initializing data
2340 2856 : this->write_nodal_data_common(fname, names, /*continuous=*/false);
2341 : #endif
2342 :
2343 1533 : if (mesh.processor_id())
2344 42 : return;
2345 :
2346 42 : int num_vars = cast_int<int>(names.size());
2347 21 : libmesh_assert_equal_to(soln.size() % num_vars, 0);
2348 273 : int num_nodes = soln.size() / num_vars;
2349 21 : libmesh_assert_equal_to(exio_helper->num_nodes, num_nodes);
2350 :
2351 : #ifndef NDEBUG
2352 21 : if (!this->get_add_sides())
2353 : {
2354 4 : int num_real_nodes = 0;
2355 4040 : for (const auto & elem : mesh.active_element_ptr_range())
2356 4036 : num_real_nodes += elem->n_nodes();
2357 4 : libmesh_assert_equal_to(num_real_nodes, num_nodes);
2358 : }
2359 : #endif
2360 :
2361 708 : for (int c=0; c<num_vars; c++)
2362 : {
2363 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2364 38 : std::vector<Real> real_parts(num_nodes);
2365 38 : std::vector<Real> imag_parts(num_nodes);
2366 : std::vector<Real> magnitudes;
2367 38 : if (_write_complex_abs)
2368 38 : magnitudes.resize(num_nodes);
2369 :
2370 : // The number of complex outputs depends on whether or not we are
2371 : // writing out the absolute values.
2372 38 : int nco = _write_complex_abs ? 3 : 2;
2373 :
2374 110896 : for (int i=0; i<num_nodes; ++i)
2375 : {
2376 110858 : real_parts[i] = soln[i*num_vars + c].real();
2377 110858 : imag_parts[i] = soln[i*num_vars + c].imag();
2378 110858 : if (_write_complex_abs)
2379 110858 : magnitudes[i] = std::abs(soln[i*num_vars + c]);
2380 : }
2381 38 : exio_helper->write_nodal_values(nco*c+1, real_parts, _timestep);
2382 38 : exio_helper->write_nodal_values(nco*c+2, imag_parts, _timestep);
2383 38 : if (_write_complex_abs)
2384 38 : exio_helper->write_nodal_values(3*c+3, magnitudes, _timestep);
2385 : #else
2386 : // Copy out this variable's solution
2387 456 : std::vector<Number> cur_soln(num_nodes);
2388 :
2389 1219856 : for (int i=0; i<num_nodes; i++)
2390 1441154 : cur_soln[i] = soln[i*num_vars + c];
2391 :
2392 418 : exio_helper->write_nodal_values(c+1,cur_soln,_timestep);
2393 : #endif
2394 : }
2395 21 : }
2396 :
2397 :
2398 :
2399 49601 : void ExodusII_IO::write_nodal_data_common(std::string fname,
2400 : const std::vector<std::string> & names,
2401 : bool continuous)
2402 : {
2403 2884 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
2404 :
2405 : // This function can be called multiple times, we only want to open
2406 : // the ExodusII file the first time it's called.
2407 49601 : if (!exio_helper->opened_for_writing)
2408 : {
2409 : // If we're appending, open() the file with read_only=false,
2410 : // otherwise create() it and write the contents of the mesh to
2411 : // it.
2412 18982 : if (_append)
2413 : {
2414 : // We do our writing only from proc 0, to avoid race
2415 : // conditions with Exodus 8
2416 365 : if (!MeshOutput<MeshBase>::mesh().processor_id())
2417 : {
2418 60 : exio_helper->open(fname.c_str(), /*read_only=*/false);
2419 : // If we're appending, it's not valid to call exio_helper->initialize()
2420 : // or exio_helper->initialize_nodal_variables(), but we do need to set up
2421 : // certain aspects of the Helper object itself, such as the number of nodes
2422 : // and elements. We do that by reading the header...
2423 60 : exio_helper->read_and_store_header_info();
2424 :
2425 : // ...and reading the block info
2426 60 : exio_helper->read_block_info();
2427 : }
2428 : // Keep other processors aware of what we've done on root
2429 : else
2430 : {
2431 295 : exio_helper->opened_for_writing = true;
2432 295 : exio_helper->current_filename = fname;
2433 : }
2434 : }
2435 : else
2436 : {
2437 36696 : exio_helper->create(fname);
2438 :
2439 : // But some of our write calls are parallel-only, due to
2440 : // calls to parallel-only getter functions.
2441 36696 : exio_helper->initialize(fname, mesh, !continuous);
2442 :
2443 18627 : exio_helper->write_nodal_coordinates(mesh, !continuous);
2444 18627 : exio_helper->write_elements(mesh, !continuous);
2445 :
2446 18627 : exio_helper->write_sidesets(mesh);
2447 18627 : exio_helper->write_nodesets(mesh);
2448 18627 : exio_helper->write_elemsets(mesh);
2449 :
2450 18627 : exio_helper->initialize_nodal_variables(names);
2451 : }
2452 : }
2453 : else
2454 : {
2455 : // We are already open for writing, so check that the filename
2456 : // passed to this function matches the filename currently in use
2457 : // by the helper.
2458 30619 : libmesh_error_msg_if(fname != exio_helper->current_filename,
2459 : "Error! This ExodusII_IO object is already associated with file: "
2460 : << exio_helper->current_filename
2461 : << ", cannot use it with requested file: "
2462 : << fname);
2463 : }
2464 49601 : }
2465 :
2466 142 : const std::vector<std::string> & ExodusII_IO::get_nodal_var_names()
2467 : {
2468 142 : exio_helper->read_var_names(ExodusII_IO_Helper::NODAL);
2469 142 : return exio_helper->nodal_var_names;
2470 : }
2471 :
2472 0 : const std::vector<std::string> & ExodusII_IO::get_elem_var_names()
2473 : {
2474 0 : exio_helper->read_var_names(ExodusII_IO_Helper::ELEMENTAL);
2475 0 : return exio_helper->elem_var_names;
2476 : }
2477 :
2478 0 : const std::vector<std::string> & ExodusII_IO::get_global_var_names()
2479 : {
2480 0 : exio_helper->read_var_names(ExodusII_IO_Helper::GLOBAL);
2481 0 : return exio_helper->global_var_names;
2482 : }
2483 :
2484 71 : const std::vector<int> & ExodusII_IO::get_elem_num_map() const
2485 : {
2486 : // We could make this function non-const and have it call
2487 : // exio_helper->read_elem_num_map() before returning a reference...
2488 : // but the intention is that this will be called some time after a
2489 : // mesh is read in, in which case it would be doing extra work to
2490 : // read in the elem_num_map twice.
2491 71 : return exio_helper->elem_num_map;
2492 : }
2493 :
2494 71 : const std::vector<int> & ExodusII_IO::get_node_num_map() const
2495 : {
2496 71 : return exio_helper->node_num_map;
2497 : }
2498 :
2499 0 : ExodusII_IO_Helper & ExodusII_IO::get_exio_helper()
2500 : {
2501 : // Provide a warning when accessing the helper object
2502 : // since it is a non-public API and is likely to see
2503 : // future API changes
2504 : libmesh_experimental();
2505 :
2506 0 : return *exio_helper;
2507 : }
2508 :
2509 :
2510 0 : void ExodusII_IO::set_hdf5_writing(bool write_hdf5)
2511 : {
2512 0 : exio_helper->set_hdf5_writing(write_hdf5);
2513 0 : }
2514 :
2515 :
2516 1988 : void ExodusII_IO::set_max_name_length(unsigned int max_length)
2517 : {
2518 1988 : exio_helper->set_max_name_length(max_length);
2519 1988 : }
2520 :
2521 :
2522 8 : void ExodusII_IO::set_discontinuous_bex(bool disc_bex)
2523 : {
2524 8 : _disc_bex = disc_bex;
2525 8 : }
2526 :
2527 :
2528 71 : void ExodusII_IO::write_discontinuous_elemental_data(bool allow_hi_order_elemental)
2529 : {
2530 71 : _allow_hi_order_elemental = allow_hi_order_elemental;
2531 71 : }
2532 :
2533 :
2534 :
2535 : // LIBMESH_HAVE_EXODUS_API is not defined, declare error() versions of functions...
2536 : #else
2537 :
2538 :
2539 :
2540 : ExodusII_IO::~ExodusII_IO () = default;
2541 :
2542 :
2543 :
2544 : void ExodusII_IO::read (const std::string &)
2545 : {
2546 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2547 : }
2548 :
2549 :
2550 :
2551 : ExodusHeaderInfo ExodusII_IO::read_header (const std::string &)
2552 : {
2553 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2554 : }
2555 :
2556 :
2557 :
2558 : void ExodusII_IO::verbose (bool)
2559 : {
2560 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2561 : }
2562 :
2563 :
2564 :
2565 : void ExodusII_IO::use_mesh_dimension_instead_of_spatial_dimension(bool)
2566 : {
2567 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2568 : }
2569 :
2570 :
2571 :
2572 : void ExodusII_IO::write_as_dimension(unsigned)
2573 : {
2574 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2575 : }
2576 :
2577 :
2578 :
2579 : void ExodusII_IO::set_coordinate_offset(Point)
2580 : {
2581 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2582 : }
2583 :
2584 :
2585 :
2586 : void ExodusII_IO::append(bool)
2587 : {
2588 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2589 : }
2590 :
2591 :
2592 :
2593 : void ExodusII_IO::write_added_sides (bool)
2594 : {
2595 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2596 : }
2597 :
2598 :
2599 :
2600 : bool ExodusII_IO::get_add_sides ()
2601 : {
2602 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2603 : return false;
2604 : }
2605 :
2606 :
2607 :
2608 : const std::vector<Real> & ExodusII_IO::get_time_steps()
2609 : {
2610 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2611 : }
2612 :
2613 :
2614 :
2615 : int ExodusII_IO::get_num_time_steps()
2616 : {
2617 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2618 : }
2619 :
2620 :
2621 : void ExodusII_IO::copy_nodal_solution(System &,
2622 : std::string,
2623 : std::string,
2624 : unsigned int)
2625 : {
2626 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2627 : }
2628 :
2629 :
2630 :
2631 : void ExodusII_IO::copy_elemental_solution(System &,
2632 : std::string,
2633 : std::string,
2634 : unsigned int)
2635 : {
2636 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2637 : }
2638 :
2639 :
2640 :
2641 : void ExodusII_IO::copy_scalar_solution(System &,
2642 : std::vector<std::string>,
2643 : std::vector<std::string>,
2644 : unsigned int)
2645 : {
2646 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2647 : }
2648 :
2649 :
2650 :
2651 : void ExodusII_IO::write_element_data (const EquationSystems &)
2652 : {
2653 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2654 : }
2655 :
2656 :
2657 :
2658 : void
2659 : ExodusII_IO::write_element_data_from_discontinuous_nodal_data
2660 : (const EquationSystems &,
2661 : const std::set<std::string> *,
2662 : const std::string & )
2663 : {
2664 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2665 : }
2666 :
2667 :
2668 :
2669 : void ExodusII_IO::write_nodal_data (const std::string &,
2670 : const std::vector<Number> &,
2671 : const std::vector<std::string> &)
2672 : {
2673 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2674 : }
2675 :
2676 :
2677 :
2678 : void ExodusII_IO::write_information_records (const std::vector<std::string> &)
2679 : {
2680 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2681 : }
2682 :
2683 :
2684 :
2685 : void ExodusII_IO::write_global_data (const std::vector<Number> &,
2686 : const std::vector<std::string> &)
2687 : {
2688 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2689 : }
2690 :
2691 :
2692 :
2693 : void ExodusII_IO::write_timestep (const std::string &,
2694 : const EquationSystems &,
2695 : const int,
2696 : const Real,
2697 : const std::set<std::string> *)
2698 : {
2699 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2700 : }
2701 :
2702 :
2703 : void ExodusII_IO::write_equation_systems (const std::string &,
2704 : const EquationSystems &,
2705 : const std::set<std::string> *)
2706 : {
2707 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2708 : }
2709 :
2710 :
2711 : void ExodusII_IO::write_elemsets()
2712 : {
2713 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2714 : }
2715 :
2716 : void
2717 : ExodusII_IO::
2718 : write_sideset_data (int,
2719 : const std::vector<std::string> &,
2720 : const std::vector<std::set<boundary_id_type>> &,
2721 : const std::vector<std::map<BoundaryInfo::BCTuple, Real>> &)
2722 : {
2723 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2724 : }
2725 :
2726 :
2727 :
2728 : void
2729 : ExodusII_IO::
2730 : read_sideset_data (int,
2731 : std::vector<std::string> &,
2732 : std::vector<std::set<boundary_id_type>> &,
2733 : std::vector<std::map<BoundaryInfo::BCTuple, Real>> &)
2734 : {
2735 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2736 : }
2737 :
2738 : void
2739 : ExodusII_IO::
2740 : get_sideset_data_indices (std::map<BoundaryInfo::BCTuple, unsigned int> &)
2741 :
2742 : {
2743 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2744 : }
2745 :
2746 : void
2747 : ExodusII_IO::
2748 : write_nodeset_data (int,
2749 : const std::vector<std::string> &,
2750 : const std::vector<std::set<boundary_id_type>> &,
2751 : const std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> &)
2752 : {
2753 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2754 : }
2755 :
2756 : void
2757 : ExodusII_IO::
2758 : read_nodeset_data (int,
2759 : std::vector<std::string> &,
2760 : std::vector<std::set<boundary_id_type>> &,
2761 : std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> &)
2762 : {
2763 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2764 : }
2765 :
2766 : void
2767 : ExodusII_IO::
2768 : get_nodeset_data_indices (std::map<BoundaryInfo::NodeBCTuple, unsigned int> &)
2769 : {
2770 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2771 : }
2772 :
2773 : void
2774 : ExodusII_IO::
2775 : write_elemset_data (int,
2776 : const std::vector<std::string> &,
2777 : const std::vector<std::set<elemset_id_type>> &,
2778 : const std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> &)
2779 : {
2780 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2781 : }
2782 :
2783 : void
2784 : ExodusII_IO::
2785 : read_elemset_data (int,
2786 : std::vector<std::string> &,
2787 : std::vector<std::set<elemset_id_type>> &,
2788 : std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> &)
2789 : {
2790 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2791 : }
2792 :
2793 : void
2794 : ExodusII_IO::
2795 : get_elemset_data_indices (std::map<std::pair<dof_id_type, elemset_id_type>, unsigned int> &)
2796 : {
2797 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2798 : }
2799 :
2800 : void ExodusII_IO::write (const std::string &)
2801 : {
2802 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2803 : }
2804 :
2805 :
2806 :
2807 : void ExodusII_IO::write_nodal_data_discontinuous (const std::string &,
2808 : const std::vector<Number> &,
2809 : const std::vector<std::string> &)
2810 : {
2811 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2812 : }
2813 :
2814 :
2815 :
2816 : void ExodusII_IO::write_nodal_data_common(std::string,
2817 : const std::vector<std::string> &,
2818 : bool)
2819 : {
2820 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2821 : }
2822 :
2823 :
2824 : const std::vector<std::string> & ExodusII_IO::get_elem_var_names()
2825 : {
2826 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2827 : }
2828 :
2829 : const std::vector<std::string> & ExodusII_IO::get_nodal_var_names()
2830 : {
2831 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2832 : }
2833 :
2834 : const std::vector<std::string> & ExodusII_IO::get_global_var_names()
2835 : {
2836 : libmesh_error_msg("ERROR, ExodusII API is not defined.");
2837 : }
2838 :
2839 : void ExodusII_IO::set_hdf5_writing(bool) {}
2840 :
2841 : #endif // LIBMESH_HAVE_EXODUS_API
2842 : } // namespace libMesh
|
