Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/IO/ExodusIIMeshReader.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief ExodusII mesh reader
9 : : \details ExodusII mesh reader class definition.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include <numeric>
14 : :
15 : : #include "NoWarning/exodusII.hpp"
16 : :
17 : : #include "ExodusIIMeshReader.hpp"
18 : : #include "ContainerUtil.hpp"
19 : : #include "Exception.hpp"
20 : : #include "UnsMesh.hpp"
21 : : #include "Reorder.hpp"
22 : :
23 : : using tk::ExodusIIMeshReader;
24 : :
25 : 22 : ExodusIIMeshReader::ExodusIIMeshReader( const std::string& filename,
26 : : int cpuwordsize,
27 : 22 : int iowordsize ) :
28 : : m_filename( filename ),
29 : : m_cpuwordsize( cpuwordsize ),
30 : : m_iowordsize( iowordsize ),
31 : : m_inFile( 0 ),
32 : : m_nnode( 0 ),
33 : : m_neblk( 0 ),
34 : : m_neset( 0 ),
35 : : m_from( 0 ),
36 : : m_till( 0 ),
37 : : m_blockid(),
38 : : m_blockid_by_type( ExoNnpe.size() ),
39 : : m_nel( ExoNnpe.size() ),
40 : : m_elemblocks(),
41 [ + - ][ + - ]: 22 : m_tri()
42 : : // *****************************************************************************
43 : : // Constructor: open Exodus II file
44 : : //! \param[in] filename File to open as ExodusII file
45 : : //! \param[in] cpuwordsize Set CPU word size, see ExodusII documentation
46 : : //! \param[in] iowordsize Set I/O word size, see ExodusII documentation
47 : : // *****************************************************************************
48 : : {
49 : : // Increase verbosity from ExodusII library in debug mode
50 : : #ifndef NDEBUG
51 [ + - ]: 22 : ex_opts( EX_DEBUG | EX_VERBOSE );
52 : : #endif
53 : :
54 : : float version;
55 : :
56 [ + - ]: 22 : m_inFile = ex_open( filename.c_str(), EX_READ, &cpuwordsize, &iowordsize,
57 : : &version );
58 : :
59 [ - + ][ - - ]: 22 : ErrChk( m_inFile > 0, "Failed to open ExodusII file: " + filename );
[ - - ][ - - ]
60 : 22 : }
61 : :
62 : 52 : ExodusIIMeshReader::~ExodusIIMeshReader() noexcept
63 : : // *****************************************************************************
64 : : // Destructor
65 : : // *****************************************************************************
66 : : {
67 [ - + ]: 52 : if ( ex_close(m_inFile) < 0 )
68 : 0 : printf( ">>> WARNING: Failed to close ExodusII file: %s\n",
69 : : m_filename.c_str() );
70 : 52 : }
71 : :
72 : : void
73 : 3 : ExodusIIMeshReader::readMesh( UnsMesh& mesh )
74 : : // *****************************************************************************
75 : : // Read ExodusII mesh file
76 : : //! \param[in] mesh Unstructured mesh object
77 : : // *****************************************************************************
78 : : {
79 : 3 : readHeader( mesh );
80 : 3 : readAllElements( mesh );
81 : 3 : readAllNodes( mesh );
82 : 3 : readSidesetFaces( mesh.bface(), mesh.faceid() );
83 : 3 : readTimeValues( mesh.vartimes() );
84 : 3 : readNodeVarNames( mesh.nodevarnames() );
85 : 6 : readNodeScalars( mesh.vartimes().size(),
86 : 3 : mesh.nodevarnames().size(),
87 : : mesh.nodevars() );
88 : 3 : }
89 : :
90 : : void
91 : 0 : ExodusIIMeshReader::readGraph( UnsMesh& mesh )
92 : : // *****************************************************************************
93 : : // Read only connectivity graph from file
94 : : //! \param[in] mesh Unstructured mesh object
95 : : // *****************************************************************************
96 : : {
97 : 0 : readHeader( mesh );
98 : 0 : readAllElements( mesh );
99 : 0 : }
100 : :
101 : : void
102 : 16 : ExodusIIMeshReader::readMeshPart(
103 : : std::vector< std::size_t >& ginpoel,
104 : : std::vector< std::size_t >& inpoel,
105 : : std::vector< std::size_t >& triinp,
106 : : std::unordered_map< std::size_t, std::size_t >& lid,
107 : : tk::UnsMesh::Coords& coord,
108 : : int numpes, int mype )
109 : : // *****************************************************************************
110 : : // Read a part of the mesh (graph and coordinates) from ExodusII file
111 : : //! \param[in,out] ginpoel Container to store element connectivity of this PE's
112 : : //! chunk of the mesh (global ids)
113 : : //! \param[in,out] inpoel Container to store element connectivity with local
114 : : //! node IDs of this PE's mesh chunk
115 : : //! \param[in,out] triinp Container to store triangle element connectivity
116 : : //! (if exists in file) with global node indices
117 : : //! \param[in,out] lid Container to store global->local node IDs of elements of
118 : : //! this PE's mesh chunk
119 : : //! \param[in,out] coord Container to store coordinates of mesh nodes of this
120 : : //! PE's mesh chunk
121 : : //! \param[in] numpes Total number of PEs (default n = 1, for a single-CPU read)
122 : : //! \param[in] mype This PE (default m = 0, for a single-CPU read)
123 : : // *****************************************************************************
124 : : {
125 [ + + ][ + - ]: 16 : Assert( mype < numpes, "Invalid input: PE id must be lower than NumPEs" );
[ + - ][ + - ]
126 [ + + ][ + + ]: 12 : Assert( ginpoel.empty() && inpoel.empty() && lid.empty() &&
[ + + ][ + + ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
127 : : coord[0].empty() && coord[1].empty() && coord[2].empty(),
128 : : "Containers to store mesh must be empty" );
129 : :
130 : : // Read info on element blocks from ExodusII file
131 [ + - ]: 8 : readElemBlockIDs();
132 : : // Get number of number of tetrahedron elements in file
133 [ + - ]: 8 : auto nel = nelem( tk::ExoElemType::TET );
134 : :
135 : : // Compute extents of element IDs of this PE's mesh chunk to read
136 : 8 : auto npes = static_cast< std::size_t >( numpes );
137 : 8 : auto pe = static_cast< std::size_t >( mype );
138 : 8 : auto chunk = nel / npes;
139 : 8 : m_from = pe * chunk;
140 : 8 : m_till = m_from + chunk;
141 [ + - ]: 8 : if (pe == npes-1) m_till += nel % npes;
142 : :
143 : : // Read tetrahedron connectivity between from and till
144 [ + - ]: 8 : readElements( {{m_from, m_till-1}}, tk::ExoElemType::TET, ginpoel );
145 : :
146 : : // Compute local data from global mesh connectivity
147 : 16 : std::vector< std::size_t > gid;
148 [ + - ]: 8 : std::tie( inpoel, gid, lid ) = tk::global2local( ginpoel );
149 : :
150 : : // Read this PE's chunk of the mesh node coordinates from file
151 [ + - ]: 8 : coord = readCoords( gid );
152 : :
153 : : // Generate set of unique faces
154 : 16 : tk::UnsMesh::FaceSet faces;
155 [ + + ]: 200 : for (std::size_t e=0; e<ginpoel.size()/4; ++e)
156 [ + + ]: 960 : for (std::size_t f=0; f<4; ++f) {
157 : 768 : const auto& tri = tk::expofa[f];
158 [ + - ]: 1536 : faces.insert( {{{ ginpoel[ e*4+tri[0] ],
159 : 768 : ginpoel[ e*4+tri[1] ],
160 : 768 : ginpoel[ e*4+tri[2] ] }}} );
161 : : }
162 : :
163 : : // Read triangle element connectivity (all triangle blocks in file)
164 [ + - ]: 8 : auto ntri = nelem( tk::ExoElemType::TRI );
165 [ + - ][ + - ]: 8 : if ( ntri !=0 ) readElements( {{0,ntri-1}}, tk::ExoElemType::TRI, triinp );
166 : :
167 : : // Keep triangles shared in (partially-read) tetrahedron mesh
168 : 16 : std::vector< std::size_t > triinp_own;
169 : 8 : std::size_t ltrid = 0; // local triangle id
170 [ + + ]: 200 : for (std::size_t e=0; e<triinp.size()/3; ++e) {
171 [ + - ]: 192 : auto i = faces.find( {{ triinp[e*3+0], triinp[e*3+1], triinp[e*3+2] }} );
172 [ + - ]: 192 : if (i != end(faces)) {
173 [ + - ]: 192 : m_tri[e] = ltrid++; // generate global->local triangle ids
174 [ + - ]: 192 : triinp_own.push_back( triinp[e*3+0] );
175 [ + - ]: 192 : triinp_own.push_back( triinp[e*3+1] );
176 [ + - ]: 192 : triinp_own.push_back( triinp[e*3+2] );
177 : : }
178 : : }
179 : 8 : triinp = std::move(triinp_own);
180 : 8 : }
181 : :
182 : : std::array< std::vector< tk::real >, 3 >
183 : 8 : ExodusIIMeshReader::readCoords( const std::vector< std::size_t >& gid ) const
184 : : // *****************************************************************************
185 : : // Read coordinates of a number of mesh nodes from ExodusII file
186 : : //! \param[in] gid Global node IDs whose coordinates to read
187 : : //! \return Vector of node coordinates read from file
188 : : // *****************************************************************************
189 : : {
190 : : // Read node coordinates from file with global node IDs given in gid
191 : 8 : return readNodes( gid );
192 : : }
193 : :
194 : : std::size_t
195 : 38 : ExodusIIMeshReader::readHeader()
196 : : // *****************************************************************************
197 : : // Read ExodusII header without setting mesh size
198 : : //! \return Number of nodes in mesh
199 : : // *****************************************************************************
200 : : {
201 : : char title[MAX_LINE_LENGTH+1];
202 : : int ndim, n, nnodeset, nelemset, nnode, neblk;
203 : :
204 [ + - ][ - + ]: 38 : ErrChk(
[ - - ][ - - ]
[ - - ]
205 : : ex_get_init( m_inFile, title, &ndim, &nnode, &n, &neblk, &nnodeset,
206 : : &nelemset ) == 0,
207 : : "Failed to read header from ExodusII file: " + m_filename );
208 : :
209 [ - + ][ - - ]: 38 : ErrChk( nnode > 0,
[ - - ][ - - ]
210 : : "Number of nodes read from ExodusII file must be larger than zero" );
211 [ - + ][ - - ]: 38 : ErrChk( neblk > 0,
[ - - ][ - - ]
212 : : "Number of element blocks read from ExodusII file must be larger "
213 : : "than zero" );
214 [ - + ][ - - ]: 38 : ErrChk( ndim == 3, "Need a 3D mesh from ExodusII file " + m_filename );
[ - - ][ - - ]
215 : :
216 : 38 : m_neblk = static_cast< std::size_t >( neblk );
217 : 38 : m_neset = static_cast< std::size_t >( nelemset );
218 : :
219 : 38 : return static_cast< std::size_t >( nnode );
220 : : }
221 : :
222 : : void
223 : 3 : ExodusIIMeshReader::readHeader( UnsMesh& mesh )
224 : : // *****************************************************************************
225 : : // Read ExodusII header with setting mesh size
226 : : //! \param[in] mesh Unstructured mesh object
227 : : // *****************************************************************************
228 : : {
229 : : // Read ExodusII file header and set mesh graph size
230 : 3 : mesh.size() = m_nnode = static_cast< std::size_t >( readHeader() );
231 : 3 : }
232 : :
233 : : void
234 : 3 : ExodusIIMeshReader::readAllNodes( UnsMesh& mesh ) const
235 : : // *****************************************************************************
236 : : // Read all node coordinates from ExodusII file
237 : : //! \param[in] mesh Unstructured mesh object
238 : : // *****************************************************************************
239 : : {
240 : 3 : mesh.x().resize( m_nnode );
241 : 3 : mesh.y().resize( m_nnode );
242 : 3 : mesh.z().resize( m_nnode );
243 : :
244 [ - + ][ - - ]: 3 : ErrChk( ex_get_coord( m_inFile, mesh.x().data(), mesh.y().data(),
[ - - ][ - - ]
245 : : mesh.z().data() ) == 0,
246 : : "Failed to read coordinates from ExodusII file: " + m_filename );
247 : 3 : }
248 : :
249 : : void
250 : 112 : ExodusIIMeshReader::readNode( std::size_t fid,
251 : : std::size_t mid,
252 : : std::vector< tk::real >& x,
253 : : std::vector< tk::real >& y,
254 : : std::vector< tk::real >& z ) const
255 : : // *****************************************************************************
256 : : // Read coordinates of a single mesh node from ExodusII file
257 : : //! \param[in] fid Node id in file whose coordinates to read
258 : : //! \param[in] mid Node id in memory to which to put new cordinates
259 : : //! \param[in,out] x Vector of x coordinates to push to
260 : : //! \param[in,out] y Vector of y coordinates to push to
261 : : //! \param[in,out] z Vector of z coordinates to push to
262 : : // *****************************************************************************
263 : : {
264 [ + - ][ + - ]: 112 : Assert( x.size() == y.size() && x.size() == z.size(), "Size mismatch" );
[ - - ][ - - ]
[ - - ]
265 [ + - ][ + - ]: 112 : Assert( mid < x.size() && mid < y.size() && mid < z.size(),
[ + - ][ - - ]
[ - - ][ - - ]
266 : : "Indexing out of bounds" );
267 : :
268 : 112 : readNode( fid, x[mid], y[mid], z[mid] );
269 : 112 : }
270 : :
271 : : void
272 : 0 : ExodusIIMeshReader::readNode( std::size_t id,
273 : : std::array< tk::real, 3 >& coord ) const
274 : : // *****************************************************************************
275 : : // Read coordinates of a single mesh node from ExodusII file
276 : : //! \param[in] id Node id whose coordinates to read
277 : : //! \param[in,out] coord Array of x, y, and z coordinates
278 : : // *****************************************************************************
279 : : {
280 : 0 : readNode( id, coord[0], coord[1], coord[2] );
281 : 0 : }
282 : :
283 : : void
284 : 112 : ExodusIIMeshReader::readNode( std::size_t id,
285 : : tk::real& x,
286 : : tk::real& y,
287 : : tk::real& z ) const
288 : : // *****************************************************************************
289 : : // Read coordinates of a single mesh node from file
290 : : //! \param[in] id Node id whose coordinates to read
291 : : //! \param[in,out] x X coordinate to write to
292 : : //! \param[in,out] y Y coordinate to write to
293 : : //! \param[in,out] z Z coordinate to write to
294 : : // *****************************************************************************
295 : : {
296 [ - + ][ - - ]: 112 : ErrChk(
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
297 : : ex_get_partial_coord( m_inFile, static_cast<int64_t>(id)+1, 1,
298 : : &x, &y, &z ) == 0,
299 : : "Failed to read coordinates of node " + std::to_string(id) +
300 : : " from ExodusII file: " + m_filename );
301 : 112 : }
302 : :
303 : : std::array< std::vector< tk::real >, 3 >
304 : 8 : ExodusIIMeshReader::readNodes( const std::vector< std::size_t >& gid ) const
305 : : // *****************************************************************************
306 : : // Read coordinates of a number of mesh nodes from ExodusII file
307 : : //! \param[in] gid Node IDs whose coordinates to read
308 : : //! \return Mesh node coordinates
309 : : // *****************************************************************************
310 : : {
311 [ + - ][ + - ]: 32 : std::vector< tk::real > px( gid.size() ), py( gid.size() ), pz( gid.size() );
[ + - ]
312 : :
313 : 8 : std::size_t i=0;
314 [ + + ][ + - ]: 120 : for (auto g : gid) readNode( g, i++, px, py, pz );
315 : :
316 : 16 : return {{ std::move(px), std::move(py), std::move(pz) }};
317 : : }
318 : :
319 : : std::size_t
320 : 28 : ExodusIIMeshReader::readElemBlockIDs()
321 : : // *****************************************************************************
322 : : // Read element block IDs from ExodusII file
323 : : //! \return Total number of nodes in mesh
324 : : // *****************************************************************************
325 : : {
326 : : // Read ExodusII file header
327 [ + - ]: 28 : auto nnode = readHeader();
328 : :
329 [ + - ]: 28 : std::vector< int > bid( m_neblk );
330 : :
331 : : // Read element block ids
332 [ + - ][ - + ]: 28 : ErrChk( ex_get_ids( m_inFile, EX_ELEM_BLOCK, bid.data()) == 0,
[ - - ][ - - ]
[ - - ]
333 : : "Failed to read element block ids from ExodusII file: " +
334 : : m_filename );
335 : :
336 : 28 : m_elemblocks.clear();
337 : 28 : m_nel.clear();
338 [ + - ]: 28 : m_nel.resize( ExoNnpe.size() );
339 : 28 : m_blockid_by_type.clear();
340 [ + - ]: 28 : m_blockid_by_type.resize( ExoNnpe.size() );
341 : :
342 : : // Fill element block ID vector
343 [ + + ]: 102 : for (auto id : bid) {
344 : : char eltype[MAX_STR_LENGTH+1];
345 : : int n, nnpe, nattr;
346 : :
347 : : // Read element block information
348 [ + - ][ - + ]: 74 : ErrChk( ex_get_block( m_inFile, EX_ELEM_BLOCK, id, eltype, &n, &nnpe,
[ - - ][ - - ]
[ - - ]
349 : : &nattr, nullptr, nullptr ) == 0,
350 : : "Failed to read element block information from ExodusII file: " +
351 : : m_filename );
352 : :
353 : : // Store ExodusII element block ID
354 [ + - ]: 74 : m_blockid.push_back( id );
355 : :
356 : 74 : auto nel = static_cast< std::size_t >( n );
357 : :
358 : : // Store info on ExodusII element blocks
359 [ + + ]: 74 : if (nnpe == 4) { // tetrahedra
360 : :
361 [ + - ]: 40 : m_elemblocks.push_back( { ExoElemType::TET, nel } );
362 : 40 : auto e = static_cast< std::size_t >( ExoElemType::TET );
363 [ + - ]: 40 : m_blockid_by_type[ e ].push_back( id );
364 [ + - ]: 40 : m_nel[ e ].push_back( nel );
365 [ - + ][ - - ]: 40 : Assert( m_blockid_by_type[e].size() == m_nel[e].size(), "Size mismatch" );
[ - - ][ - - ]
366 : :
367 [ + + ]: 34 : } else if (nnpe == 3) { // triangles
368 : :
369 [ + - ]: 23 : m_elemblocks.push_back( { ExoElemType::TRI, nel } );
370 : 23 : auto e = static_cast< std::size_t >( ExoElemType::TRI );
371 [ + - ]: 23 : m_blockid_by_type[ e ].push_back( id );
372 [ + - ]: 23 : m_nel[ e ].push_back( nel );
373 [ - + ][ - - ]: 23 : Assert( m_blockid_by_type[e].size() == m_nel[e].size(), "Size mismatch" );
[ - - ][ - - ]
374 : :
375 : : }
376 : : }
377 : :
378 : 56 : return nnode;
379 : : }
380 : :
381 : : void
382 : 3 : ExodusIIMeshReader::readAllElements( UnsMesh& mesh )
383 : : // *****************************************************************************
384 : : // Read all element blocks and mesh connectivity from ExodusII file
385 : : //! \param[inout] mesh Unstructured mesh object to store mesh in
386 : : // *****************************************************************************
387 : : {
388 : : // Read element block ids
389 : 3 : readElemBlockIDs();
390 : :
391 [ + + ]: 12 : for (auto id : m_blockid) {
392 : : char eltype[MAX_STR_LENGTH+1];
393 : : int nel, nnpe, nattr;
394 : :
395 : : // Read element block information
396 [ + - ][ - + ]: 9 : ErrChk( ex_get_block( m_inFile, EX_ELEM_BLOCK, id, eltype, &nel, &nnpe,
[ - - ][ - - ]
[ - - ]
397 : : &nattr, nullptr, nullptr ) == 0,
398 : : "Failed to read element block information from ExodusII file: " +
399 : : m_filename );
400 : :
401 : : // Read element connectivity
402 : 9 : auto connectsize = static_cast< std::size_t >( nel*nnpe );
403 [ + + ]: 9 : if (nnpe == 4) { // tetrahedra
404 : :
405 [ + - ]: 14 : std::vector< int > inpoel( connectsize );
406 [ + - ][ - + ]: 7 : ErrChk( ex_get_conn( m_inFile, EX_ELEM_BLOCK, id, inpoel.data(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
407 : : nullptr, nullptr ) == 0,
408 : : "Failed to read " + std::string(eltype) + " element connectivity from "
409 : : "ExodusII file: " + m_filename );
410 [ + + ]: 4299 : for (auto n : inpoel)
411 [ + - ]: 4292 : mesh.tetinpoel().push_back( static_cast< std::size_t >( n ) );
412 : :
413 [ + + ]: 2 : } else if (nnpe == 3) { // triangles
414 : :
415 [ + - ]: 2 : std::vector< int > inpoel( connectsize );
416 [ + - ][ - + ]: 1 : ErrChk( ex_get_conn( m_inFile, EX_ELEM_BLOCK, id, inpoel.data(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
417 : : nullptr, nullptr ) == 0,
418 : : "Failed to read " + std::string(eltype) + " element connectivity from "
419 : : "ExodusII file: " + m_filename );
420 [ + + ]: 73 : for (auto n : inpoel)
421 [ + - ]: 72 : mesh.triinpoel().push_back( static_cast< std::size_t >( n ) );
422 : :
423 : : }
424 : : }
425 : :
426 : : // Shift node IDs to start from zero
427 : 3 : shiftToZero( mesh.triinpoel() );
428 : 3 : shiftToZero( mesh.tetinpoel() );
429 : 3 : }
430 : :
431 : : void
432 : 5656 : ExodusIIMeshReader::readElements( const std::array< std::size_t, 2 >& ext,
433 : : tk::ExoElemType elemtype,
434 : : std::vector< std::size_t >& conn ) const
435 : : // *****************************************************************************
436 : : // Read element connectivity of a number of mesh cells from ExodusII file
437 : : //! \param[in] ext Extents of element IDs whose connectivity to read:
438 : : //! [from...till), using zero-based element IDs, where 'from' >=0, inclusive
439 : : //! and 'till < 'maxelements', where 'maxelements' is the total number of
440 : : //! elements of all element blocks in the file of the requested cell type.
441 : : //! Note that 'maxelements' can be queried by nelem().
442 : : //! \param[in] elemtype Element type
443 : : //! \param[inout] conn Connectivity vector to push to
444 : : //! \note Must be preceded by a call to readElemBlockIDs()
445 : : //! \details This function takes the extents of element IDs in a zero-based
446 : : //! fashion. These input extents can be thought of "absolute" extents that
447 : : //! denote lowest and the largest-1 element IDs to be read from file.
448 : : // *****************************************************************************
449 : : {
450 [ - + ][ - - ]: 5656 : Assert( tk::sumsize(m_blockid_by_type) > 0,
[ - - ][ - - ]
451 : : "A call to this function must be preceded by a call to "
452 : : "ExodusIIMeshReader::readElemBlockIDs()" );
453 [ + - ][ + - ]: 5656 : Assert( ext[0] <= ext[1] &&
[ + - ][ + - ]
[ + - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
454 : : ext[0] < nelem(elemtype) &&
455 : : ext[1] < nelem(elemtype),
456 : : "Invalid element ID extents. Of the requested extents [from...till), "
457 : : "'from' must be lower than or equal to 'till', and they must be in "
458 : : "the range [0...maxelements), where 'maxelements' is the total "
459 : : "number of elements of all element blocks in the file of the "
460 : : "requested cell type. Requested element ID extents: ["
461 : : + std::to_string(ext[0]) + "..." + std::to_string(ext[1])
462 : : + "), 'maxelements' of cell type with "
463 : : + std::to_string( ExoNnpe[ static_cast<std::size_t>(elemtype) ] )
464 : : + " nodes per cell in file '" + m_filename + "': "
465 : : + std::to_string( nelem( elemtype ) ) );
466 : :
467 : 5656 : auto e = static_cast< std::size_t >( elemtype );
468 : : // List of number of elements of all blocks of element type requested
469 : 5656 : const auto& nel = m_nel[e];
470 : : // List of element block IDs for element type requested
471 : 5656 : const auto& bid = m_blockid_by_type[e];
472 : :
473 : : // Compute lower and upper element block ids to read from based on extents
474 : 5656 : std::size_t lo_bid = 0, hi_bid = 0, offset = 0;
475 [ + + ]: 32736 : for (std::size_t b=0; b<nel.size(); ++b) {
476 : 27080 : std::size_t lo = offset; // lo (min) elem ID in block
477 : 27080 : std::size_t hi = offset + nel[b] - 1; // hi (max) elem ID in block
478 [ + + ][ + + ]: 27080 : if (ext[0] >= lo && ext[0] <= hi) lo_bid = b;
[ + + ]
479 [ + + ][ + + ]: 27080 : if (ext[1] >= lo && ext[1] <= hi) hi_bid = b;
[ + + ]
480 : 27080 : offset += nel[b];
481 : : }
482 : :
483 [ + - ][ + - ]: 5656 : Assert( lo_bid < nel.size() && lo_bid < bid.size(),
[ - - ][ - - ]
[ - - ]
484 : : "Invalid start block ID" );
485 [ + - ][ + - ]: 5656 : Assert( hi_bid < nel.size() && hi_bid < bid.size(),
[ - - ][ - - ]
[ - - ]
486 : : "Invalid end block ID" );
487 : :
488 : : // Compute relative extents based on absolute ones for each block to read from
489 : 11312 : std::vector< std::array< std::size_t, 2 > > rext;
490 : 5656 : offset = 0;
491 [ + + ]: 12104 : for (std::size_t b=0; b<lo_bid; ++b) offset += nel[b];
492 [ + + ]: 19840 : for (std::size_t b=lo_bid; b<=hi_bid; ++b) {
493 : 14184 : std::size_t lo = offset;
494 : 14184 : std::size_t hi = offset + nel[b] - 1;
495 : 14184 : std::size_t le = 1, he = nel[b];
496 [ + + ][ + - ]: 14184 : if (ext[0] >= lo && ext[0] <= hi) le = ext[0] - lo + 1;
[ + + ]
497 [ + - ][ + + ]: 14184 : if (ext[1] >= lo && ext[1] <= hi) he = ext[1] - lo + 1;
[ + + ]
498 [ + - ][ + - ]: 14184 : Assert( le >= 1 && le <= nel[b] && he >= 1 && he <= nel[b],
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ]
499 : : "Relative index out of block" );
500 [ + - ]: 14184 : rext.push_back( {{ le, he }} );
501 : 14184 : offset += nel[b];
502 : : }
503 : :
504 [ + - ][ - + ]: 14184 : Assert( std::accumulate(
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
505 : : std::next(rext.cbegin()), rext.cend(), rext[0][1]-rext[0][0]+1,
506 : : []( std::size_t n, const std::array< std::size_t, 2 >& r )
507 : : { return n + r[1] - r[0] + 1; }
508 : : ) == ext[1]-ext[0]+1,
509 : : "Total number of elements to read incorrect, requested extents: " +
510 : : std::to_string(ext[0]) + " ... " + std::to_string(ext[1]) );
511 : :
512 : 11312 : std::vector< int > inpoel;
513 : :
514 : : // Read element connectivity from file
515 : 5656 : std::size_t B = 0;
516 [ + + ]: 19840 : for (auto b=lo_bid; b<=hi_bid; ++b, ++B) {
517 : 14184 : const auto& r = rext[B];
518 [ + - ]: 28368 : std::vector< int > c( (r[1]-r[0]+1) * ExoNnpe[e] );
519 [ + - ][ - + ]: 14184 : ErrChk( ex_get_partial_conn( m_inFile,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
520 : : EX_ELEM_BLOCK,
521 : : bid[b],
522 : : static_cast< int64_t >( r[0] ),
523 : : static_cast< int64_t >( r[1]-r[0]+1 ),
524 : : c.data(),
525 : : nullptr,
526 : : nullptr ) == 0,
527 : : "Failed to read element connectivity of elements [" +
528 : : std::to_string(r[0]) + "..." + std::to_string(r[1]) +
529 : : "] from element block " + std::to_string(bid[b]) + " in ExodusII "
530 : : "file: " + m_filename );
531 [ + - ]: 14184 : inpoel.reserve( inpoel.size() + c.size() );
532 [ + - ][ + - ]: 14184 : std::move( begin(c), end(c), std::back_inserter(inpoel) );
533 : : }
534 : :
535 [ - + ][ - - ]: 5656 : Assert( inpoel.size() == (ext[1]-ext[0]+1)*ExoNnpe[e],
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
536 : : "Failed to read element connectivity of elements [" +
537 : : std::to_string(ext[0]) + "..." + std::to_string(ext[1]) + ") from "
538 : : "ExodusII file: " + m_filename );
539 : :
540 : : // Put in element connectivity using zero-based node indexing
541 [ + + ]: 7394742 : for (auto& i : inpoel) --i;
542 [ + - ]: 5656 : conn.reserve( conn.size() + inpoel.size() );
543 [ + - ][ + - ]: 5656 : std::move( begin(inpoel), end(inpoel), std::back_inserter(conn) );
544 : 5656 : }
545 : :
546 : : void
547 : 8 : ExodusIIMeshReader::readFaces( std::vector< std::size_t >& conn ) const
548 : : // *****************************************************************************
549 : : // Read face connectivity of a number of boundary faces from ExodusII file
550 : : //! \param[inout] conn Connectivity vector to push to
551 : : //! \details This function reads in the total number of boundary faces,
552 : : //! also called triangle-elements in the EXO2 file, and their connectivity.
553 : : // *****************************************************************************
554 : : {
555 : : // Return quietly if no triangle elements in file
556 [ - + ]: 8 : if (nelem(tk::ExoElemType::TRI) == 0) return;
557 : :
558 : : // Read triangle boundary-face connectivity (all the triangle element block)
559 [ + - ]: 8 : readElements( {{0,nelem(tk::ExoElemType::TRI)-1}}, tk::ExoElemType::TRI,
560 : : conn );
561 : : }
562 : :
563 : : std::vector< std::size_t >
564 : 1 : ExodusIIMeshReader::readNodemap()
565 : : // *****************************************************************************
566 : : // Read local to global node-ID map from ExodusII file
567 : : //! \return node_map Vector mapping the local Exodus node-IDs to global node-IDs
568 : : //! \details The node-map is required to get the "Exodus-global" node-IDs from
569 : : //! the "Exodus-internal" node-IDs, which are returned from the exodus APIs.
570 : : //! The node-IDs in the exodus file are referred to as the "Exodus-global"
571 : : //! node-IDs or "fileIDs" in Quinoa.
572 : : // *****************************************************************************
573 : : {
574 : : // Read triangle boundary-face connectivity
575 [ + - ]: 1 : auto nnode = readElemBlockIDs();
576 : :
577 : : // Create array to store node-number map
578 [ + - ]: 2 : std::vector< int > node_map( nnode );
579 : :
580 : : // Read in the node number map to map the above nodes to the global node-IDs
581 [ + - ][ - + ]: 1 : ErrChk( ex_get_id_map( m_inFile, EX_NODE_MAP, node_map.data() ) == 0,
[ - - ][ - - ]
[ - - ]
582 : : "Failed to read node map length from ExodusII file: " );
583 : :
584 [ + - ]: 1 : std::vector< std::size_t > node_map1( nnode );
585 : :
586 [ + + ]: 15 : for (std::size_t i=0; i<nnode; ++i)
587 : : {
588 : 14 : node_map1[i] = static_cast< std::size_t >(node_map[i]-1);
589 : : }
590 : :
591 : 2 : return node_map1;
592 : : }
593 : :
594 : : std::map< int, std::vector< std::size_t > >
595 : 7 : ExodusIIMeshReader::readSidesetNodes()
596 : : // *****************************************************************************
597 : : // Read node list of all side sets from ExodusII file
598 : : //! \return Node lists mapped to side set ids
599 : : // *****************************************************************************
600 : : {
601 : : // Read ExodusII file header (fills m_neset)
602 : 7 : readHeader();
603 : :
604 : : // Node lists mapped to side set ids
605 : 7 : std::map< int, std::vector< std::size_t > > side;
606 : :
607 [ + - ]: 7 : if (m_neset > 0) {
608 : : // Read all side set ids from file
609 [ + - ]: 14 : std::vector< int > ids( m_neset );
610 [ + - ][ - + ]: 7 : ErrChk( ex_get_ids( m_inFile, EX_SIDE_SET, ids.data() ) == 0,
[ - - ][ - - ]
[ - - ]
611 : : "Failed to read side set ids from ExodusII file: " + m_filename );
612 : : // Read in node list for all side sets
613 [ + + ]: 28 : for (auto i : ids) {
614 : : int nface, nnode;
615 : : // Read number of faces and number of distribution factors in side set i
616 [ + - ][ - + ]: 21 : ErrChk( ex_get_set_param( m_inFile, EX_SIDE_SET, i, &nface, &nnode ) == 0,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
617 : : "Failed to read side set " + std::to_string(i) + " parameters "
618 : : "from ExodusII file: " + m_filename );
619 : : // Read number of nodes in side set i (overwrite nnode)
620 [ + - ][ - + ]: 21 : ErrChk( ex_get_side_set_node_list_len( m_inFile, i, &nnode ) == 0,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
621 : : "Failed to read side set " + std::to_string(i) + " node list "
622 : : "length from ExodusII file: " + m_filename );
623 [ - + ][ - - ]: 21 : Assert(nnode > 0, "Number of nodes = 0 in side set" + std::to_string(i));
[ - - ][ - - ]
[ - - ]
624 [ + - ]: 42 : std::vector< int > df( static_cast< std::size_t >( nface ) );
625 [ + - ]: 42 : std::vector< int > nodes( static_cast< std::size_t >( nnode ) );
626 : : // Read in node list for side set i
627 [ + - ][ - + ]: 21 : ErrChk( ex_get_side_set_node_list( m_inFile, i, df.data(), nodes.data() )
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
628 : : == 0, "Failed to read node list of side set " +
629 : : std::to_string(i) + " from ExodusII file: " +
630 : : m_filename );
631 : : // Make node list unique
632 [ + - ]: 21 : tk::unique( nodes );
633 : : // Store 0-based node ID list as std::size_t vector instead of ints
634 [ + - ]: 21 : auto& list = side[ i ];
635 [ + + ][ + - ]: 9576 : for (auto n : nodes) list.push_back( static_cast<std::size_t>(n-1) );
636 : : }
637 : : }
638 : :
639 : 7 : return side;
640 : : }
641 : :
642 : : void
643 : 14 : ExodusIIMeshReader::readSidesetFaces(
644 : : std::map< int, std::vector< std::size_t > >& bface,
645 : : std::map< int, std::vector< std::size_t > >& faces )
646 : : // *****************************************************************************
647 : : // Read side sets from ExodusII file
648 : : //! \param[in,out] bface Elem ids of side sets to read into
649 : : //! \param[in,out] faces Elem-relative face ids of tets of side sets
650 : : // *****************************************************************************
651 : : {
652 : : // Read element block ids
653 : 14 : readElemBlockIDs();
654 : :
655 [ + + ]: 14 : if (m_neset > 0) {
656 : : // Read side set ids from file
657 [ + - ]: 22 : std::vector< int > ids( m_neset );
658 [ + - ][ - + ]: 11 : ErrChk( ex_get_ids( m_inFile, EX_SIDE_SET, ids.data() ) == 0,
[ - - ][ - - ]
[ - - ]
659 : : "Failed to read side set ids from ExodusII file: " + m_filename );
660 : :
661 : : // Read all side sets from file
662 [ + + ]: 36 : for (auto i : ids) {
663 : : int nface, nnode;
664 : :
665 : : // Read number of faces in side set
666 [ + - ][ - + ]: 25 : ErrChk( ex_get_set_param( m_inFile, EX_SIDE_SET, i, &nface, &nnode ) == 0,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
667 : : "Failed to read side set " + std::to_string(i) + " parameters "
668 : : "from ExodusII file: " + m_filename );
669 : :
670 [ - + ][ - - ]: 25 : Assert(nface > 0, "Number of faces = 0 in side set" + std::to_string(i));
[ - - ][ - - ]
[ - - ]
671 : :
672 [ + - ]: 50 : std::vector< int > exoelem( static_cast< std::size_t >( nface ) );
673 [ + - ]: 50 : std::vector< int > exoface( static_cast< std::size_t >( nface ) );
674 : :
675 : : // Read in file-internal element ids and relative face ids for side set
676 [ + - ][ - + ]: 25 : ErrChk( ex_get_set( m_inFile, EX_SIDE_SET, i, exoelem.data(),
[ - - ][ - - ]
[ - - ][ - - ]
677 : : exoface.data() ) == 0,
678 : : "Failed to read side set " + std::to_string(i) );
679 : :
680 : : // Store file-internal element ids of side set
681 [ + - ]: 25 : auto& elem = bface[i];
682 [ + - ]: 25 : elem.resize( exoelem.size() );
683 : 25 : std::size_t j = 0;
684 [ + + ]: 16907 : for (auto e : exoelem) elem[j++] = static_cast< std::size_t >( e-1 );
685 : :
686 : : // Store zero-based relative face ids of side set
687 [ + - ]: 25 : auto& face = faces[i];
688 [ + - ]: 25 : face.resize( exoface.size() );
689 : 25 : j = 0;
690 [ + + ]: 16907 : for (auto n : exoface) face[j++] = static_cast< std::size_t >( n-1 );
691 : :
692 [ + - ][ - + ]: 16907 : Assert( std::all_of( begin(face), end(face),
[ - - ][ - - ]
[ - - ]
693 : : [](std::size_t f){ return f<4; } ),
694 : : "Relative face id of side set must be between 0 and 3" );
695 [ - + ][ - - ]: 25 : Assert( elem.size() == face.size(), "Size mismatch" );
[ - - ][ - - ]
696 : : }
697 : : }
698 : 14 : }
699 : :
700 : : std::pair< tk::ExoElemType, std::size_t >
701 : 0 : ExodusIIMeshReader::blkRelElemId( std::size_t id ) const
702 : : // *****************************************************************************
703 : : // Compute element-block-relative element id and element type
704 : : //! \param[in] id (ExodusII) file-internal element id
705 : : //! \return Element type the internal id points to and element id relative to
706 : : //! cell-type
707 : : //! \details This function takes an internal element id, which in general can
708 : : //! point to any element block in the ExodusII file and thus we do not know
709 : : //! which element type a block contains. It then computes which cell type the
710 : : //! id points to and computes the relative index for the given cell type. This
711 : : //! is necessary because elements are read in from file by from potentially
712 : : //! multiple blocks by cell type.
713 : : //! \note Must be preceded by a call to readElemBlockIDs()
714 : : // *****************************************************************************
715 : : {
716 : 0 : auto TRI = tk::ExoElemType::TRI;
717 : 0 : auto TET = tk::ExoElemType::TET;
718 : :
719 : 0 : std::size_t e = 0; // counts elements (independent of cell type)
720 : 0 : std::size_t ntri = 0; // counts triangle elements
721 : 0 : std::size_t ntet = 0; // counts tetrahedron elements
722 : :
723 [ - - ]: 0 : for (const auto& b : m_elemblocks) { // walk all element blocks in order
724 : 0 : e += b.second; // increment file-internal element id
725 [ - - ]: 0 : if (e > id) { // found element block for internal id
726 [ - - ]: 0 : if (b.first == TRI) { // if triangle block
727 : 0 : return { TRI, id-ntet }; // return cell type and triangle id
728 [ - - ]: 0 : } else if (b.first == TET) { // if tetrahedron block
729 : 0 : return { TET, id-ntri }; // return cell type and tetrahedron id
730 : : }
731 : : }
732 : : // increment triangle and tetrahedron elements independently
733 [ - - ]: 0 : if (b.first == TRI)
734 : 0 : ntri += b.second;
735 [ - - ]: 0 : else if (b.first == TET)
736 : 0 : ntet += b.second;
737 : : }
738 : :
739 [ - - ][ - - ]: 0 : Throw( " Exodus internal element id not found" );
[ - - ]
740 : : }
741 : :
742 : : std::vector< std::size_t >
743 : 0 : ExodusIIMeshReader::triinpoel(
744 : : std::map< int, std::vector< std::size_t > >& belem,
745 : : const std::map< int, std::vector< std::size_t > >& faces,
746 : : const std::vector< std::size_t >& ginpoel,
747 : : const std::vector< std::size_t >& triinp ) const
748 : : // *****************************************************************************
749 : : // Generate triangle face connectivity for side sets
750 : : //! \param[in,out] belem File-internal elem ids of side sets
751 : : //! \param[in] faces Elem-relative face ids of side sets
752 : : //! \param[in] ginpoel Tetrahedron element connectivity with global nodes
753 : : //! \param[in] triinp Triangle element connectivity with global nodes
754 : : //! (if exists in file)
755 : : //! \return Triangle face connectivity with global node IDs of side sets
756 : : //! \details This function takes lists of file-internal element ids (in belem)
757 : : //! for side sets and does two things: (1) generates face connectivity (with
758 : : //! global node IDs) for side sets, and (2) converts the (ExodusII)
759 : : //! file-internal element IDs to face ids so that they can be used to index
760 : : //! into the face connectivity. The IDs in belem are modified and the face
761 : : //! connectivity (for boundary faces only) is returned.
762 : : //! \note Must be preceded by a call to readElemBlockIDs()
763 : : // *****************************************************************************
764 : : {
765 [ - - ][ - - ]: 0 : Assert( !(m_from == 0 && m_till == 0),
[ - - ][ - - ]
[ - - ]
766 : : "Lower and upper tetrahedron id bounds must not both be zero" );
767 : :
768 : : // This will contain one of our final results: face (triangle) connectivity
769 : : // for the side sets only. The difference between bnd_triinpoel and triinpoel
770 : : // is that triinpoel is a triangle element connectivity, independent of side
771 : : // sets, while bnd_triinpoel is a triangle connectivity only for side sets.
772 : 0 : std::vector< std::size_t > bnd_triinpoel;
773 : :
774 : : // Storage for boundary face lists for each side set on this PE
775 : 0 : std::map< int, std::vector< std::size_t > > belem_own;
776 : :
777 : 0 : std::size_t f = 0; // counts all faces
778 [ - - ]: 0 : for (auto& ss : belem) { // for all side sets
779 : :
780 : : // insert side set id into new map
781 [ - - ]: 0 : auto& b = belem_own[ ss.first ];
782 : : // get element-relative face ids for side set
783 [ - - ]: 0 : const auto& face = tk::cref_find( faces, ss.first );
784 : 0 : std::size_t s = 0; // counts side set faces
785 [ - - ]: 0 : for (auto& i : ss.second) { // for all faces on side set
786 : :
787 : : // compute element-block-relative element id and element type
788 [ - - ]: 0 : auto r = blkRelElemId( i );
789 : :
790 : : // extract boundary face connectivity based on element type
791 : 0 : bool localface = false;
792 [ - - ]: 0 : if (r.first == tk::ExoElemType::TRI) {
793 : :
794 [ - - ]: 0 : auto t = m_tri.find(r.second);
795 [ - - ]: 0 : if (t != end(m_tri)) { // only if triangle id exists on this PE
796 [ - - ][ - - ]: 0 : Assert( t->second < triinp.size()/3,
[ - - ][ - - ]
797 : : "Indexing out of triangle connectivity" );
798 : : // generate triangle (face) connectivity using global node ids
799 [ - - ]: 0 : bnd_triinpoel.push_back( triinp[ t->second*3 + 0 ] );
800 [ - - ]: 0 : bnd_triinpoel.push_back( triinp[ t->second*3 + 1 ] );
801 [ - - ]: 0 : bnd_triinpoel.push_back( triinp[ t->second*3 + 2 ] );
802 : 0 : localface = true;
803 : : }
804 : :
805 [ - - ]: 0 : } else if (r.first == tk::ExoElemType::TET) {
806 : :
807 [ - - ][ - - ]: 0 : if (r.second >= m_from && r.second < m_till) { // if tet is on this PE
808 : 0 : auto t = r.second - m_from;
809 [ - - ][ - - ]: 0 : Assert( t < ginpoel.size()/4,
[ - - ][ - - ]
810 : : "Indexing out of tetrahedron connectivity" );
811 : : // get ExodusII face-node numbering for side sets, see ExodusII
812 : : // manual figure on "Sideset side Numbering"
813 : 0 : const auto& tri = tk::expofa[ face[s] ];
814 : : // generate triangle (face) connectivity using global node ids, note
815 : : // the switched node order, 0,2,1, as lpofa is different from expofa
816 [ - - ]: 0 : bnd_triinpoel.push_back( ginpoel[ t*4 + tri[0] ] );
817 [ - - ]: 0 : bnd_triinpoel.push_back( ginpoel[ t*4 + tri[1] ] );
818 [ - - ]: 0 : bnd_triinpoel.push_back( ginpoel[ t*4 + tri[2] ] );
819 : 0 : localface = true;
820 : : }
821 : :
822 : : }
823 : :
824 : 0 : ++s;
825 : :
826 : : // generate PE-local face id for side set (this is to be used to index
827 : : // into triinpoel)
828 [ - - ][ - - ]: 0 : if (localface) b.push_back( f++ );
829 : : }
830 : :
831 : : // if no faces on this side set (on this PE), remove side set id
832 [ - - ][ - - ]: 0 : if (b.empty()) belem_own.erase( ss.first );
833 : : }
834 : :
835 : 0 : belem = std::move(belem_own);
836 : :
837 : 0 : return bnd_triinpoel;
838 : : }
839 : :
840 : : void
841 : 3 : ExodusIIMeshReader::readNodeVarNames( std::vector< std::string >& nv ) const
842 : : // *****************************************************************************
843 : : // Read the names of nodal output variables from ExodusII file
844 : : //! \param[in,out] nv Nodal variable names
845 : : // *****************************************************************************
846 : : {
847 : : #if defined(__clang__)
848 : : #pragma clang diagnostic push
849 : : #pragma clang diagnostic ignored "-Wvla"
850 : : #pragma clang diagnostic ignored "-Wvla-extension"
851 : : #elif defined(STRICT_GNUC)
852 : : #pragma GCC diagnostic push
853 : : #pragma GCC diagnostic ignored "-Wvla"
854 : : #endif
855 : :
856 : 3 : int numvars = 0;
857 : :
858 [ + - ][ - + ]: 3 : ErrChk(
[ - - ][ - - ]
[ - - ]
859 : : ex_get_variable_param( m_inFile, EX_NODE_BLOCK, &numvars ) == 0,
860 : : "Failed to read nodal output variable parameters from ExodusII file: " +
861 : : m_filename );
862 : :
863 [ - + ]: 3 : if (numvars) {
864 : :
865 : 0 : char* names[ static_cast< std::size_t >( numvars ) ];
866 [ - - ]: 0 : for (int i=0; i<numvars; ++i)
867 : 0 : names[i] = static_cast<char*>( calloc((MAX_STR_LENGTH+1), sizeof(char)) );
868 : :
869 [ - - ][ - - ]: 0 : ErrChk( ex_get_variable_names( m_inFile,
[ - - ][ - - ]
[ - - ]
870 : : EX_NODAL,
871 : : numvars,
872 : : names ) == 0,
873 : : "Failed to read nodal variable names from ExodusII file: " +
874 : : m_filename );
875 : :
876 [ - - ]: 0 : nv.resize( static_cast< std::size_t >( numvars ) );
877 : 0 : std::size_t i = 0;
878 [ - - ][ - - ]: 0 : for (auto& n : nv) n = names[ i++ ];
879 : :
880 : : }
881 : :
882 : : #if defined(__clang__)
883 : : #pragma clang diagnostic pop
884 : : #elif defined(STRICT_GNUC)
885 : : #pragma GCC diagnostic pop
886 : : #endif
887 : 3 : }
888 : :
889 : : void
890 : 3 : ExodusIIMeshReader::readTimeValues( std::vector< tk::real >& tv ) const
891 : : // *****************************************************************************
892 : : // Read time values from ExodusII file
893 : : //! \param[in] tv Vector of time values at which field data is saved
894 : : // *****************************************************************************
895 : : {
896 : : auto num_time_steps =
897 : 3 : static_cast< std::size_t >( ex_inquire_int( m_inFile, EX_INQ_TIME ) );
898 : :
899 [ - + ]: 3 : if (num_time_steps) {
900 [ - - ]: 0 : tv.resize( num_time_steps, 0.0 );
901 [ - - ][ - - ]: 0 : ErrChk( ex_get_all_times( m_inFile, tv.data() ) == 0,
[ - - ][ - - ]
902 : : "Failed to read time values from ExodusII file: " + m_filename );
903 : : }
904 : 3 : }
905 : :
906 : : void
907 : 3 : ExodusIIMeshReader::readNodeScalars(
908 : : std::size_t ntime,
909 : : std::size_t nvar,
910 : : std::vector< std::vector< std::vector< tk::real > > >& var ) const
911 : : // *****************************************************************************
912 : : // Read node scalar fields from ExodusII file
913 : : //! \param[in] ntime Number of time steps to read
914 : : //! \param[in] nvar Number of variables to read
915 : : //! \param[in] var Vector of nodal variables to read to: inner vector: nodes,
916 : : //! middle vector: (physics) variable, outer vector: time step
917 : : // *****************************************************************************
918 : : {
919 : 3 : var.resize( ntime );
920 [ - + ]: 3 : for (auto& v : var) {
921 [ - - ]: 0 : v.resize( nvar );
922 [ - - ][ - - ]: 0 : for (auto& n : v) n.resize( m_nnode );
923 : : }
924 : :
925 [ - + ]: 3 : for (std::size_t t=0; t<var.size(); ++t) {
926 [ - - ]: 0 : for (std::size_t id=0; id<var[t].size(); ++id) {
927 [ - - ][ - - ]: 0 : ErrChk( ex_get_var( m_inFile,
[ - - ][ - - ]
928 : : static_cast< int >( t+1 ),
929 : : EX_NODAL,
930 : : static_cast< int >( id+1 ),
931 : : 1,
932 : : static_cast< int64_t >( var[t][id].size() ),
933 : : var[t][id].data() ) == 0,
934 : : "Failed to read node scalar from ExodusII file: " + m_filename );
935 : : }
936 : : }
937 : 3 : }
938 : :
939 : : std::size_t
940 : 11347 : ExodusIIMeshReader::nelem( tk::ExoElemType elemtype ) const
941 : : // *****************************************************************************
942 : : // Return number of elements in all mesh blocks for a given elem type in file
943 : : //! \param[in] elemtype Element type
944 : : //! \return Number of elements in all blocks for the elem type
945 : : //! \note Must be preceded by a call to readElemBlockIDs()
946 : : // *****************************************************************************
947 : : {
948 : 11347 : auto e = static_cast< std::size_t >( elemtype );
949 : 11347 : return std::accumulate( m_nel[e].cbegin(), m_nel[e].cend(), 0u );
950 : : }
|
