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