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