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