Quinoa regression test code coverage report
Current view: top level - IO - ExodusIIMeshReader.cpp (source / functions) Hit Total Coverage
Commit: Quinoa_v0.3-957-gb4f0efae0 Lines: 291 312 93.3 %
Date: 2021-11-11 13:17:06 Functions: 23 28 82.1 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 262 838 31.3 %

           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                 :        927 : ExodusIIMeshReader::ExodusIIMeshReader( const std::string& filename,
      26                 :            :                                         int cpuwordsize,
      27                 :        927 :                                         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 [ +  - ][ +  - ]:        927 :   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         [ +  - ]:        927 :   ex_opts( EX_DEBUG | EX_VERBOSE );
      52                 :            :   #endif
      53                 :            : 
      54                 :            :   float version;
      55                 :            : 
      56         [ +  - ]:        927 :   m_inFile = ex_open( filename.c_str(), EX_READ, &cpuwordsize, &iowordsize,
      57                 :            :                       &version );
      58                 :            : 
      59 [ -  + ][ -  - ]:        927 :   ErrChk( m_inFile > 0, "Failed to open ExodusII file: " + filename );
         [ -  - ][ -  - ]
      60                 :        927 : }
      61                 :            : 
      62                 :       2715 : ExodusIIMeshReader::~ExodusIIMeshReader() noexcept
      63                 :            : // *****************************************************************************
      64                 :            : //  Destructor
      65                 :            : // *****************************************************************************
      66                 :            : {
      67         [ -  + ]:       2715 :   if ( ex_close(m_inFile) < 0 )
      68                 :          0 :     printf( ">>> WARNING: Failed to close ExodusII file: %s\n",
      69                 :            :             m_filename.c_str() );
      70                 :       2715 : }
      71                 :            : 
      72                 :            : void
      73                 :         33 : ExodusIIMeshReader::readMesh( UnsMesh& mesh )
      74                 :            : // *****************************************************************************
      75                 :            : //  Read ExodusII mesh file
      76                 :            : //! \param[in] mesh Unstructured mesh object
      77                 :            : // *****************************************************************************
      78                 :            : {
      79                 :         33 :   readHeader( mesh );
      80                 :         33 :   readAllElements( mesh );
      81                 :         33 :   readAllNodes( mesh );
      82                 :         33 :   readSidesetFaces( mesh.bface(), mesh.faceid() );
      83                 :         33 :   readTimeValues( mesh.vartimes() );
      84                 :         33 :   readNodeVarNames( mesh.nodevarnames() );
      85                 :         66 :   readNodeScalars( mesh.vartimes().size(),
      86                 :         33 :                    mesh.nodevarnames().size(),
      87                 :            :                    mesh.nodevars() );
      88                 :         33 : }
      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                 :        663 : 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 [ -  + ][ -  - ]:        663 :   Assert( mype < numpes, "Invalid input: PE id must be lower than NumPEs" );
         [ -  - ][ -  - ]
     126 [ +  - ][ +  - ]:        663 :   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         [ +  - ]:        663 :   readElemBlockIDs();
     132                 :            :   // Get number of number of tetrahedron elements in file
     133         [ +  - ]:        663 :   auto nel = nelem( tk::ExoElemType::TET );
     134                 :            : 
     135                 :            :   // Compute extents of element IDs of this PE's mesh chunk to read
     136                 :        663 :   auto npes = static_cast< std::size_t >( numpes );
     137                 :        663 :   auto pe = static_cast< std::size_t >( mype );
     138                 :        663 :   auto chunk = nel / npes;
     139                 :        663 :   m_from = pe * chunk;
     140                 :        663 :   m_till = m_from + chunk;
     141         [ +  + ]:        663 :   if (pe == npes-1) m_till += nel % npes;
     142                 :            : 
     143                 :            :   // Read tetrahedron connectivity between from and till
     144         [ +  - ]:        663 :   readElements( {{m_from, m_till-1}}, tk::ExoElemType::TET, ginpoel );
     145                 :            : 
     146                 :            :   // Compute local data from global mesh connectivity
     147                 :       1326 :   std::vector< std::size_t > gid;
     148         [ +  - ]:        663 :   std::tie( inpoel, gid, lid ) = tk::global2local( ginpoel );
     149                 :            : 
     150                 :            :   // Read this PE's chunk of the mesh node coordinates from file
     151         [ +  - ]:        663 :   coord = readCoords( gid );
     152                 :            : 
     153                 :            :   // Generate set of unique faces
     154                 :       1326 :   tk::UnsMesh::FaceSet faces;
     155         [ +  + ]:    2303339 :   for (std::size_t e=0; e<ginpoel.size()/4; ++e)
     156         [ +  + ]:   11513380 :     for (std::size_t f=0; f<4; ++f) {
     157                 :    9210704 :       const auto& tri = tk::expofa[f];
     158         [ +  - ]:   18421408 :       faces.insert( {{{ ginpoel[ e*4+tri[0] ],
     159                 :    9210704 :                         ginpoel[ e*4+tri[1] ],
     160                 :    9210704 :                         ginpoel[ e*4+tri[2] ] }}} );
     161                 :            :     }
     162                 :            : 
     163                 :            :   // Read triangle element connectivity (all triangle blocks in file)
     164         [ +  - ]:        663 :   auto ntri = nelem( tk::ExoElemType::TRI );
     165 [ +  + ][ +  - ]:        663 :   if ( ntri !=0 ) readElements( {{0,ntri-1}}, tk::ExoElemType::TRI, triinp );
     166                 :            : 
     167                 :            :   // Keep triangles shared in (partially-read) tetrahedron mesh
     168                 :       1326 :   std::vector< std::size_t > triinp_own;
     169                 :        663 :   std::size_t ltrid = 0;        // local triangle id
     170         [ +  + ]:    1893537 :   for (std::size_t e=0; e<triinp.size()/3; ++e) {
     171         [ +  - ]:    1892874 :     auto i = faces.find( {{ triinp[e*3+0], triinp[e*3+1], triinp[e*3+2] }} );
     172         [ +  + ]:    1892874 :     if (i != end(faces)) {
     173         [ +  - ]:     558750 :       m_tri[e] = ltrid++;       // generate global->local triangle ids
     174         [ +  - ]:     558750 :       triinp_own.push_back( triinp[e*3+0] );
     175         [ +  - ]:     558750 :       triinp_own.push_back( triinp[e*3+1] );
     176         [ +  - ]:     558750 :       triinp_own.push_back( triinp[e*3+2] );
     177                 :            :     }
     178                 :            :   }
     179                 :        663 :   triinp = std::move(triinp_own);
     180                 :        663 : }
     181                 :            : 
     182                 :            : std::array< std::vector< tk::real >, 3 >
     183                 :        663 : 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                 :        663 :   return readNodes( gid );
     192                 :            : }
     193                 :            : 
     194                 :            : std::size_t
     195                 :       1359 : 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 [ +  - ][ -  + ]:       1359 :   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 [ -  + ][ -  - ]:       1359 :   ErrChk( nnode > 0,
         [ -  - ][ -  - ]
     210                 :            :           "Number of nodes read from ExodusII file must be larger than zero" );
     211 [ -  + ][ -  - ]:       1359 :   ErrChk( neblk > 0,
         [ -  - ][ -  - ]
     212                 :            :           "Number of element blocks read from ExodusII file must be larger "
     213                 :            :           "than zero" );
     214 [ -  + ][ -  - ]:       1359 :   ErrChk( ndim == 3, "Need a 3D mesh from ExodusII file " + m_filename );
         [ -  - ][ -  - ]
     215                 :            : 
     216                 :       1359 :   m_neblk = static_cast< std::size_t >( neblk );
     217                 :       1359 :   m_neset = static_cast< std::size_t >( nelemset );
     218                 :            : 
     219                 :       1359 :   return static_cast< std::size_t >( nnode );
     220                 :            : }
     221                 :            : 
     222                 :            : void
     223                 :         33 : 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                 :         33 :   mesh.size() = m_nnode = static_cast< std::size_t >( readHeader() );
     231                 :         33 : }
     232                 :            : 
     233                 :            : void
     234                 :         33 : ExodusIIMeshReader::readAllNodes( UnsMesh& mesh ) const
     235                 :            : // *****************************************************************************
     236                 :            : //  Read all node coordinates from ExodusII file
     237                 :            : //! \param[in] mesh Unstructured mesh object
     238                 :            : // *****************************************************************************
     239                 :            : {
     240                 :         33 :   mesh.x().resize( m_nnode );
     241                 :         33 :   mesh.y().resize( m_nnode );
     242                 :         33 :   mesh.z().resize( m_nnode );
     243                 :            : 
     244 [ -  + ][ -  - ]:         33 :   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                 :         33 : }
     248                 :            : 
     249                 :            : void
     250                 :    1440005 : 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 [ +  - ][ +  - ]:    1440005 :   Assert( x.size() == y.size() && x.size() == z.size(), "Size mismatch" );
         [ -  - ][ -  - ]
                 [ -  - ]
     265 [ +  - ][ +  - ]:    1440005 :   Assert( mid < x.size() && mid < y.size() && mid < z.size(),
         [ +  - ][ -  - ]
         [ -  - ][ -  - ]
     266                 :            :           "Indexing out of bounds" );
     267                 :            : 
     268                 :    1440005 :   readNode( fid, x[mid], y[mid], z[mid] );
     269                 :    1440005 : }
     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                 :    1440005 : 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 [ -  + ][ -  - ]:    1440005 :   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                 :    1440005 : }
     302                 :            : 
     303                 :            : std::array< std::vector< tk::real >, 3 >
     304                 :        663 : 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 [ +  - ][ +  - ]:       2652 :   std::vector< tk::real > px( gid.size() ), py( gid.size() ), pz( gid.size() );
                 [ +  - ]
     312                 :            : 
     313                 :        663 :   std::size_t i=0;
     314 [ +  + ][ +  - ]:    1440668 :   for (auto g : gid) readNode( g, i++, px, py, pz );
     315                 :            : 
     316                 :       1326 :   return {{ std::move(px), std::move(py), std::move(pz) }};
     317                 :            : }
     318                 :            : 
     319                 :            : std::size_t
     320                 :        960 : 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         [ +  - ]:        960 :   auto nnode = readHeader();
     328                 :            : 
     329         [ +  - ]:        960 :   std::vector< int > bid( m_neblk );
     330                 :            : 
     331                 :            :   // Read element block ids
     332 [ +  - ][ -  + ]:        960 :   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                 :        960 :   m_elemblocks.clear();
     337                 :        960 :   m_nel.clear();
     338         [ +  - ]:        960 :   m_nel.resize( ExoNnpe.size() );
     339                 :        960 :   m_blockid_by_type.clear();
     340         [ +  - ]:        960 :   m_blockid_by_type.resize( ExoNnpe.size() );
     341                 :            : 
     342                 :            :   // Fill element block ID vector
     343         [ +  + ]:       2781 :   for (auto id : bid) {
     344                 :            :     char eltype[MAX_STR_LENGTH+1];
     345                 :            :     int n, nnpe, nattr;
     346                 :            : 
     347                 :            :     // Read element block information
     348 [ +  - ][ -  + ]:       1821 :     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         [ +  - ]:       1821 :     m_blockid.push_back( id );
     355                 :            : 
     356                 :       1821 :     auto nel = static_cast< std::size_t >( n );
     357                 :            : 
     358                 :            :     // Store info on ExodusII element blocks
     359         [ +  + ]:       1821 :     if (nnpe == 4) {        // tetrahedra
     360                 :            : 
     361         [ +  - ]:        908 :       m_elemblocks.push_back( { ExoElemType::TET, nel } );
     362                 :        908 :       auto e = static_cast< std::size_t >( ExoElemType::TET );
     363         [ +  - ]:        908 :       m_blockid_by_type[ e ].push_back( id );
     364         [ +  - ]:        908 :       m_nel[ e ].push_back( nel );
     365 [ -  + ][ -  - ]:        908 :       Assert( m_blockid_by_type[e].size() == m_nel[e].size(), "Size mismatch" );
         [ -  - ][ -  - ]
     366                 :            : 
     367         [ +  + ]:        913 :     } else if (nnpe == 3) { // triangles
     368                 :            : 
     369         [ +  - ]:        909 :       m_elemblocks.push_back( { ExoElemType::TRI, nel } );
     370                 :        909 :       auto e = static_cast< std::size_t >( ExoElemType::TRI );
     371         [ +  - ]:        909 :       m_blockid_by_type[ e ].push_back( id );
     372         [ +  - ]:        909 :       m_nel[ e ].push_back( nel );
     373 [ -  + ][ -  - ]:        909 :       Assert( m_blockid_by_type[e].size() == m_nel[e].size(), "Size mismatch" );
         [ -  - ][ -  - ]
     374                 :            : 
     375                 :            :     }
     376                 :            :   }
     377                 :            : 
     378                 :       1920 :   return nnode;
     379                 :            : }
     380                 :            : 
     381                 :            : void
     382                 :         33 : 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                 :         33 :   readElemBlockIDs();
     390                 :            : 
     391         [ +  + ]:         74 :   for (auto id : m_blockid) {
     392                 :            :     char eltype[MAX_STR_LENGTH+1];
     393                 :            :     int nel, nnpe, nattr;
     394                 :            : 
     395                 :            :     // Read element block information
     396 [ +  - ][ -  + ]:         41 :     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                 :         41 :     auto connectsize = static_cast< std::size_t >( nel*nnpe );
     403         [ +  + ]:         41 :     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         [ +  + ]:     196479 :       for (auto n : inpoel)
     411         [ +  - ]:     196472 :         mesh.tetinpoel().push_back( static_cast< std::size_t >( n ) );
     412                 :            : 
     413         [ +  + ]:         34 :     } else if (nnpe == 3) {    // triangles
     414                 :            : 
     415         [ +  - ]:         64 :       std::vector< int > inpoel( connectsize );
     416 [ +  - ][ -  + ]:         32 :       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         [ +  + ]:     191291 :       for (auto n : inpoel)
     421         [ +  - ]:     191259 :         mesh.triinpoel().push_back( static_cast< std::size_t >( n ) );
     422                 :            : 
     423                 :            :     }
     424                 :            :   }
     425                 :            : 
     426                 :            :   // Shift node IDs to start from zero
     427                 :         33 :   shiftToZero( mesh.triinpoel() );
     428                 :         33 :   shiftToZero( mesh.tetinpoel() );
     429                 :         33 : }
     430                 :            : 
     431                 :            : void
     432                 :       1288 : 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 [ -  + ][ -  - ]:       1288 :   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 [ +  - ][ +  - ]:       1288 :   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                 :       1288 :   auto e = static_cast< std::size_t >( elemtype );
     468                 :            :   // List of number of elements of all blocks of element type requested
     469                 :       1288 :   const auto& nel = m_nel[e];
     470                 :            :   // List of element block IDs for element type requested
     471                 :       1288 :   const auto& bid = m_blockid_by_type[e];
     472                 :            : 
     473                 :            :   // Compute lower and upper element block ids to read from based on extents
     474                 :       1288 :   std::size_t lo_bid = 0, hi_bid = 0, offset = 0;
     475         [ +  + ]:       2576 :   for (std::size_t b=0; b<nel.size(); ++b) {
     476                 :       1288 :     std::size_t lo = offset;                    // lo (min) elem ID in block
     477                 :       1288 :     std::size_t hi = offset + nel[b] - 1;       // hi (max) elem ID in block
     478 [ +  - ][ +  - ]:       1288 :     if (ext[0] >= lo && ext[0] <= hi) lo_bid = b;
                 [ +  - ]
     479 [ +  - ][ +  - ]:       1288 :     if (ext[1] >= lo && ext[1] <= hi) hi_bid = b;
                 [ +  - ]
     480                 :       1288 :     offset += nel[b];
     481                 :            :   }
     482                 :            : 
     483 [ +  - ][ +  - ]:       1288 :   Assert( lo_bid < nel.size() && lo_bid < bid.size(),
         [ -  - ][ -  - ]
                 [ -  - ]
     484                 :            :           "Invalid start block ID" );
     485 [ +  - ][ +  - ]:       1288 :   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                 :       2576 :   std::vector< std::array< std::size_t, 2 > > rext;
     490                 :       1288 :   offset = 0;
     491         [ -  + ]:       1288 :   for (std::size_t b=0; b<lo_bid; ++b) offset += nel[b];
     492         [ +  + ]:       2576 :   for (std::size_t b=lo_bid; b<=hi_bid; ++b) {
     493                 :       1288 :     std::size_t lo = offset;
     494                 :       1288 :     std::size_t hi = offset + nel[b] - 1;
     495                 :       1288 :     std::size_t le = 1, he = nel[b];
     496 [ +  - ][ +  - ]:       1288 :     if (ext[0] >= lo && ext[0] <= hi) le = ext[0] - lo + 1;
                 [ +  - ]
     497 [ +  - ][ +  - ]:       1288 :     if (ext[1] >= lo && ext[1] <= hi) he = ext[1] - lo + 1;
                 [ +  - ]
     498 [ +  - ][ +  - ]:       1288 :     Assert( le >= 1 && le <= nel[b] && he >= 1 && he <= nel[b],
         [ +  - ][ +  - ]
         [ -  - ][ -  - ]
                 [ -  - ]
     499                 :            :             "Relative index out of block" );
     500         [ +  - ]:       1288 :     rext.push_back( {{ le, he }} );
     501                 :       1288 :     offset += nel[b];
     502                 :            :   }
     503                 :            : 
     504 [ +  - ][ -  + ]:       1288 :   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                 :       2576 :   std::vector< int > inpoel;
     513                 :            : 
     514                 :            :   // Read element connectivity from file
     515                 :       1288 :   std::size_t B = 0;
     516         [ +  + ]:       2576 :   for (auto b=lo_bid; b<=hi_bid; ++b, ++B) {
     517                 :       1288 :     const auto& r = rext[B];
     518         [ +  - ]:       2576 :     std::vector< int > c( (r[1]-r[0]+1) * ExoNnpe[e] );
     519 [ +  - ][ -  + ]:       1288 :     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         [ +  - ]:       1288 :     inpoel.reserve( inpoel.size() + c.size() );
     532 [ +  - ][ +  - ]:       1288 :     std::move( begin(c), end(c), std::back_inserter(inpoel) );
     533                 :            :   }
     534                 :            : 
     535 [ -  + ][ -  - ]:       1288 :   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         [ +  + ]:   14890614 :   for (auto& i : inpoel) --i;
     542         [ +  - ]:       1288 :   conn.reserve( conn.size() + inpoel.size() );
     543 [ +  - ][ +  - ]:       1288 :   std::move( begin(inpoel), end(inpoel), std::back_inserter(conn) );
     544                 :       1288 : }
     545                 :            : 
     546                 :            : void
     547                 :          0 : 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         [ -  - ]:          0 :   if (nelem(tk::ExoElemType::TRI) == 0) return;
     557                 :            : 
     558                 :            :   // Read triangle boundary-face connectivity (all the triangle element block)
     559         [ -  - ]:          0 :   readElements( {{0,nelem(tk::ExoElemType::TRI)-1}}, tk::ExoElemType::TRI,
     560                 :            :                 conn );
     561                 :            : }
     562                 :            : 
     563                 :            : std::vector< std::size_t >
     564                 :          0 : 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         [ -  - ]:          0 :   auto nnode = readElemBlockIDs();
     576                 :            : 
     577                 :            :   // Create array to store node-number map
     578         [ -  - ]:          0 :   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 [ -  - ][ -  - ]:          0 :   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         [ -  - ]:          0 :   std::vector< std::size_t > node_map1( nnode );
     585                 :            : 
     586         [ -  - ]:          0 :   for (std::size_t i=0; i<nnode; ++i)
     587                 :            :   {
     588                 :          0 :           node_map1[i] = static_cast< std::size_t >(node_map[i]-1);
     589                 :            :   }
     590                 :            : 
     591                 :          0 :   return node_map1;
     592                 :            : }
     593                 :            : 
     594                 :            : std::map< int, std::vector< std::size_t > >
     595                 :        135 : 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                 :        135 :   readHeader();
     603                 :            : 
     604                 :            :   // Node lists mapped to side set ids
     605                 :        135 :   std::map< int, std::vector< std::size_t > > side;
     606                 :            : 
     607         [ +  + ]:        135 :   if (m_neset > 0) {
     608                 :            :     // Read all side set ids from file
     609         [ +  - ]:        214 :     std::vector< int > ids( m_neset );
     610 [ +  - ][ -  + ]:        107 :     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         [ +  + ]:        549 :     for (auto i : ids) {
     614                 :            :       int nface, nnode;
     615                 :            :       // Read number of faces and number of distribution factors in side set i
     616 [ +  - ][ -  + ]:        442 :       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 [ +  - ][ -  + ]:        442 :       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 [ -  + ][ -  - ]:        442 :       Assert(nnode > 0, "Number of nodes = 0 in side set" + std::to_string(i));
         [ -  - ][ -  - ]
                 [ -  - ]
     624         [ +  - ]:        884 :       std::vector< int > df( static_cast< std::size_t >( nface ) );
     625         [ +  - ]:        884 :       std::vector< int > nodes( static_cast< std::size_t >( nnode ) );
     626                 :            :       // Read in node list for side set i
     627 [ +  - ][ -  + ]:        442 :       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         [ +  - ]:        442 :       tk::unique( nodes );
     633                 :            :       // Store 0-based node ID list as std::size_t vector instead of ints
     634         [ +  - ]:        442 :       auto& list = side[ i ];
     635 [ +  + ][ +  - ]:     167411 :       for (auto n : nodes) list.push_back( static_cast<std::size_t>(n-1) );
     636                 :            :     }
     637                 :            :   }
     638                 :            : 
     639                 :        135 :   return side;
     640                 :            : }
     641                 :            : 
     642                 :            : void
     643                 :        264 : 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                 :        264 :   readElemBlockIDs();
     654                 :            : 
     655         [ +  + ]:        264 :   if (m_neset > 0) {
     656                 :            :     // Read side set ids from file
     657         [ +  - ]:        406 :     std::vector< int > ids( m_neset );
     658 [ +  - ][ -  + ]:        203 :     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         [ +  + ]:        995 :     for (auto i : ids) {
     663                 :            :       int nface, nnode;
     664                 :            : 
     665                 :            :       // Read number of faces in side set
     666 [ +  - ][ -  + ]:        792 :       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 [ -  + ][ -  - ]:        792 :       Assert(nface > 0, "Number of faces = 0 in side set" + std::to_string(i));
         [ -  - ][ -  - ]
                 [ -  - ]
     671                 :            : 
     672         [ +  - ]:       1584 :       std::vector< int > exoelem( static_cast< std::size_t >( nface ) );
     673         [ +  - ]:       1584 :       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 [ +  - ][ -  + ]:        792 :       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         [ +  - ]:        792 :       auto& elem = bface[i];
     682         [ +  - ]:        792 :       elem.resize( exoelem.size() );
     683                 :        792 :       std::size_t j = 0;
     684         [ +  + ]:     533642 :       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         [ +  - ]:        792 :       auto& face = faces[i];
     688         [ +  - ]:        792 :       face.resize( exoface.size() );
     689                 :        792 :       j = 0;
     690         [ +  + ]:     533642 :       for (auto n : exoface) face[j++] = static_cast< std::size_t >( n-1 );
     691                 :            : 
     692 [ +  - ][ -  + ]:     533642 :       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 [ -  + ][ -  - ]:        792 :       Assert( elem.size() == face.size(), "Size mismatch" );
         [ -  - ][ -  - ]
     696                 :            :     }
     697                 :            :   }
     698                 :        264 : }
     699                 :            : 
     700                 :            : std::pair< tk::ExoElemType, std::size_t >
     701                 :     968524 : 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                 :     968524 :   auto TRI = tk::ExoElemType::TRI;
     717                 :     968524 :   auto TET = tk::ExoElemType::TET;
     718                 :            : 
     719                 :     968524 :   std::size_t e = 0;            // counts elements (independent of cell type)
     720                 :     968524 :   std::size_t ntri = 0;         // counts triangle elements
     721                 :     968524 :   std::size_t ntet = 0;         // counts tetrahedron elements
     722                 :            : 
     723         [ +  - ]:     998784 :   for (const auto& b : m_elemblocks) {  // walk all element blocks in order
     724                 :     998784 :     e += b.second;                      // increment file-internal element id
     725         [ +  + ]:     998784 :     if (e > id) {                       // found element block for internal id
     726         [ +  + ]:     968524 :       if (b.first == TRI) {             // if triangle block
     727                 :     904792 :         return { TRI, id-ntet };        // return cell type and triangle id
     728         [ +  - ]:      63732 :       } else if (b.first == TET) {      // if tetrahedron block
     729                 :      63732 :         return { TET, id-ntri };        // return cell type and tetrahedron id
     730                 :            :       }
     731                 :            :     }
     732                 :            :     // increment triangle and tetrahedron elements independently
     733         [ -  + ]:      30260 :     if (b.first == TRI)
     734                 :          0 :       ntri += b.second;
     735         [ +  - ]:      30260 :     else if (b.first == TET)
     736                 :      30260 :       ntet += b.second;
     737                 :            :   }
     738                 :            : 
     739 [ -  - ][ -  - ]:          0 :   Throw( " Exodus internal element id not found" );
                 [ -  - ]
     740                 :            : }
     741                 :            : 
     742                 :            : std::vector< std::size_t >
     743                 :        663 : 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 [ +  + ][ -  + ]:        663 :   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                 :        663 :   std::vector< std::size_t > bnd_triinpoel;
     773                 :            : 
     774                 :            :   // Storage for boundary face lists for each side set on this PE
     775                 :       1326 :   std::map< int, std::vector< std::size_t > > belem_own;
     776                 :            : 
     777                 :        663 :   std::size_t f = 0;            // counts all faces
     778         [ +  + ]:       2616 :   for (auto& ss : belem) {      // for all side sets
     779                 :            : 
     780                 :            :     // insert side set id into new map
     781         [ +  - ]:       1953 :     auto& b = belem_own[ ss.first ];
     782                 :            :     // get element-relative face ids for side set
     783         [ +  - ]:       1953 :     const auto& face = tk::cref_find( faces, ss.first );
     784                 :       1953 :     std::size_t s = 0;          // counts side set faces
     785         [ +  + ]:     970477 :     for (auto& i : ss.second) { // for all faces on side set
     786                 :            : 
     787                 :            :       // compute element-block-relative element id and element type
     788         [ +  - ]:     968524 :       auto r = blkRelElemId( i );
     789                 :            : 
     790                 :            :       // extract boundary face connectivity based on element type
     791                 :     968524 :       bool localface = false;
     792         [ +  + ]:     968524 :       if (r.first == tk::ExoElemType::TRI) {
     793                 :            : 
     794         [ +  - ]:     904792 :         auto t = m_tri.find(r.second);
     795         [ +  + ]:     904792 :         if (t != end(m_tri)) {  // only if triangle id exists on this PE
     796 [ -  + ][ -  - ]:     291236 :           Assert( t->second < triinp.size()/3,
         [ -  - ][ -  - ]
     797                 :            :                   "Indexing out of triangle connectivity" );
     798                 :            :           // generate triangle (face) connectivity using global node ids
     799         [ +  - ]:     291236 :           bnd_triinpoel.push_back( triinp[ t->second*3 + 0 ] );
     800         [ +  - ]:     291236 :           bnd_triinpoel.push_back( triinp[ t->second*3 + 1 ] );
     801         [ +  - ]:     291236 :           bnd_triinpoel.push_back( triinp[ t->second*3 + 2 ] );
     802                 :     291236 :           localface = true;
     803                 :            :         }
     804                 :            : 
     805         [ +  - ]:      63732 :       } else if (r.first == tk::ExoElemType::TET) {
     806                 :            : 
     807 [ +  + ][ +  + ]:      63732 :         if (r.second >= m_from && r.second < m_till) {  // if tet is on this PE
     808                 :      16038 :           auto t = r.second - m_from;
     809 [ -  + ][ -  - ]:      16038 :           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                 :      16038 :           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         [ +  - ]:      16038 :           bnd_triinpoel.push_back( ginpoel[ t*4 + tri[0] ] );
     817         [ +  - ]:      16038 :           bnd_triinpoel.push_back( ginpoel[ t*4 + tri[1] ] );
     818         [ +  - ]:      16038 :           bnd_triinpoel.push_back( ginpoel[ t*4 + tri[2] ] );
     819                 :      16038 :           localface = true;
     820                 :            :         }
     821                 :            : 
     822                 :            :       }
     823                 :            : 
     824                 :     968524 :       ++s;
     825                 :            : 
     826                 :            :       // generate PE-local face id for side set (this is to be used to index
     827                 :            :       // into triinpoel)
     828 [ +  + ][ +  - ]:     968524 :       if (localface) b.push_back( f++ );
     829                 :            :     }
     830                 :            : 
     831                 :            :     // if no faces on this side set (on this PE), remove side set id
     832 [ +  + ][ +  - ]:       1953 :     if (b.empty()) belem_own.erase( ss.first );
     833                 :            :   }
     834                 :            : 
     835                 :        663 :   belem = std::move(belem_own);
     836                 :            : 
     837                 :       1326 :   return bnd_triinpoel;
     838                 :            : }
     839                 :            : 
     840                 :            : void
     841                 :         33 : 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                 :         33 :   int numvars = 0;
     857                 :            : 
     858 [ +  - ][ -  + ]:         33 :   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         [ +  + ]:         33 :   if (numvars) {
     864                 :            : 
     865                 :         30 :     char* names[ static_cast< std::size_t >( numvars ) ];
     866         [ +  + ]:        210 :     for (int i=0; i<numvars; ++i)
     867                 :        180 :       names[i] = static_cast<char*>( calloc((MAX_STR_LENGTH+1), sizeof(char)) );
     868                 :            : 
     869 [ +  - ][ -  + ]:         30 :     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         [ +  - ]:         30 :     nv.resize( static_cast< std::size_t >( numvars ) );
     877                 :         30 :     std::size_t i = 0;
     878 [ +  + ][ +  - ]:        240 :     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                 :         33 : }
     888                 :            : 
     889                 :            : void
     890                 :         33 : 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                 :         33 :     static_cast< std::size_t >( ex_inquire_int( m_inFile, EX_INQ_TIME ) );
     898                 :            : 
     899         [ +  + ]:         33 :   if (num_time_steps) {
     900         [ +  - ]:         30 :     tv.resize( num_time_steps, 0.0 );
     901 [ -  + ][ -  - ]:         30 :     ErrChk( ex_get_all_times( m_inFile, tv.data() ) == 0,
         [ -  - ][ -  - ]
     902                 :            :              "Failed to read time values from ExodusII file: " + m_filename );
     903                 :            :   }
     904                 :         33 : }
     905                 :            : 
     906                 :            : void
     907                 :         33 : 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                 :         33 :   var.resize( ntime );
     920         [ +  + ]:        123 :   for (auto& v : var) {
     921         [ +  - ]:         90 :     v.resize( nvar );
     922 [ +  + ][ +  - ]:        630 :     for (auto& n : v) n.resize( m_nnode );
     923                 :            :   }
     924                 :            : 
     925         [ +  + ]:        123 :   for (std::size_t t=0; t<var.size(); ++t) {
     926         [ +  + ]:        630 :     for (std::size_t id=0; id<var[t].size(); ++id) {
     927 [ -  + ][ -  - ]:        540 :       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                 :         33 : }
     938                 :            : 
     939                 :            : std::size_t
     940                 :       3902 : 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                 :       3902 :   auto e = static_cast< std::size_t >( elemtype );
     949                 :       3902 :   return std::accumulate( m_nel[e].cbegin(), m_nel[e].cend(), 0u );
     950                 :            : }

Generated by: LCOV version 1.14