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