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