Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/IO/MeshFactory.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     Unstructured mesh reader and writer factory
       9                 :            :   \details   Unstructured mesh reader and writer factory.
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : 
      13                 :            : #include <string>
      14                 :            : 
      15                 :            : #include "MeshFactory.hpp"
      16                 :            : #include "MeshDetect.hpp"
      17                 :            : #include "Timer.hpp"
      18                 :            : #include "Reader.hpp"
      19                 :            : #include "GmshMeshReader.hpp"
      20                 :            : #include "NetgenMeshReader.hpp"
      21                 :            : #include "ExodusIIMeshReader.hpp"
      22                 :            : #include "HyperMeshReader.hpp"
      23                 :            : #include "UGRIDMeshReader.hpp"
      24                 :            : #include "RDGFLOMeshReader.hpp"
      25                 :            : #include "ASCMeshReader.hpp"
      26                 :            : #include "NetgenMeshWriter.hpp"
      27                 :            : #include "GmshMeshWriter.hpp"
      28                 :            : #include "ExodusIIMeshWriter.hpp"
      29                 :            : #include "DerivedData.hpp"
      30                 :            : #include "Reorder.hpp"
      31                 :            : #include "QuinoaBuildConfig.hpp"
      32                 :            : 
      33                 :            : #ifdef HAS_OMEGA_H
      34                 :            :   #include "Omega_h_MeshReader.hpp"
      35                 :            : #endif
      36                 :            : 
      37                 :            : namespace tk {
      38                 :            : 
      39                 :            : UnsMesh
      40                 :         47 : readUnsMesh( const tk::Print& print,
      41                 :            :              const std::string& filename,
      42                 :            :              std::pair< std::string, tk::real >& timestamp )
      43                 :            : // *****************************************************************************
      44                 :            : //  Read unstructured mesh from file
      45                 :            : //! \param[in] print Pretty printer
      46                 :            : //! \param[in] filename Filename to read mesh from
      47                 :            : //! \param[out] timestamp A time stamp consisting of a timer label (a string),
      48                 :            : //!   and a time state (a tk::real in seconds) measuring the mesh read time
      49                 :            : //! \return Unstructured mesh object
      50                 :            : // *****************************************************************************
      51                 :            : {
      52 [ +  - ][ +  - ]:         47 :   print.diagstart( "Reading mesh from file '" + filename + "' ..." );
                 [ +  - ]
      53                 :            : 
      54                 :            :   // Read in mesh
      55                 :         47 :   tk::Timer t;
      56                 :            :  
      57                 :            :   //! Create unstructured mesh to store mesh
      58                 :         47 :   UnsMesh mesh;
      59                 :            : 
      60         [ +  - ]:         47 :   const auto meshtype = detectInput( filename );
      61                 :            : 
      62         [ +  + ]:         47 :   if (meshtype == MeshReaderType::GMSH)
      63 [ +  - ][ +  - ]:          5 :     GmshMeshReader( filename ).readMesh( mesh );
      64         [ +  + ]:         42 :   else if (meshtype == MeshReaderType::NETGEN)
      65 [ +  - ][ +  - ]:          2 :     NetgenMeshReader( filename ).readMesh( mesh );
      66         [ +  + ]:         40 :   else if (meshtype == MeshReaderType::EXODUSII)
      67 [ +  - ][ +  - ]:         33 :     ExodusIIMeshReader( filename ).readMesh( mesh );
      68         [ +  + ]:          7 :   else if (meshtype == MeshReaderType::ASC)
      69 [ +  - ][ +  - ]:          4 :     ASCMeshReader( filename ).readMesh( mesh );
      70         [ +  + ]:          3 :   else if (meshtype == MeshReaderType::UGRID)
      71 [ +  - ][ +  - ]:          1 :     UGRIDMeshReader( filename ).readMesh( mesh );
      72         [ +  + ]:          2 :   else if (meshtype == MeshReaderType::RDGFLO)
      73 [ +  - ][ +  - ]:          1 :     RDGFLOMeshReader( filename ).readMesh( mesh );
      74         [ +  - ]:          1 :   else if (meshtype == MeshReaderType::HYPER)
      75 [ +  - ][ +  - ]:          1 :     HyperMeshReader( filename ).readMesh( mesh );
      76                 :            : 
      77                 :            :   timestamp =
      78 [ +  - ][ +  - ]:         47 :     std::make_pair( "Read mesh from file '" + filename + '\'', t.dsec() );
         [ +  - ][ +  - ]
      79                 :            : 
      80 [ +  - ][ +  - ]:         47 :   print.diagend( "done" );
      81                 :            : 
      82                 :            :   // Return (move out) mesh object
      83                 :         94 :   return mesh;
      84                 :            : }
      85                 :            : 
      86                 :            : std::vector< std::pair< std::string, tk::real > >
      87                 :         18 : writeUnsMesh( const tk::Print& print,
      88                 :            :               const std::string& filename,
      89                 :            :               UnsMesh& mesh,
      90                 :            :               bool reorder )
      91                 :            : // *****************************************************************************
      92                 :            : //  Write unstructured mesh to file
      93                 :            : //! \param[in] print Pretty printer
      94                 :            : //! \param[in] filename Filename to write mesh to
      95                 :            : //! \param[in] mesh Unstructured mesh object to write from
      96                 :            : //! \param[in] reorder Whether to also reorder mesh nodes
      97                 :            : //! \return Vector of time stamps consisting of a timer label (a string), and a
      98                 :            : //!   time state (a tk::real in seconds) measuring the renumber and the mesh
      99                 :            : //!   write time
     100                 :            : // *****************************************************************************
     101                 :            : {
     102                 :         18 :   std::vector< std::pair< std::string, tk::real > > times;
     103                 :            : 
     104                 :         18 :   tk::Timer t;
     105                 :            : 
     106                 :            :   // If mesh has tetrahedra but no triangles, generate triangle connectivity
     107 [ +  + ][ +  + ]:         18 :   if (!mesh.tetinpoel().empty() && mesh.triinpoel().empty()) {
                 [ +  + ]
     108 [ +  - ][ +  - ]:          6 :     print.diagstart( "Generating missing surface mesh ..." );
     109                 :            : 
     110                 :          6 :     const auto& inpoel = mesh.tetinpoel();        // get tet connectivity
     111         [ +  - ]:         12 :     auto esup = tk::genEsup( inpoel, 4 );         // elements surrounding points
     112         [ +  - ]:         12 :     auto esuel = tk::genEsuelTet( inpoel, esup ); // elems surrounding elements
     113                 :          6 :     auto& triinpoel = mesh.triinpoel();
     114                 :            :     // collect boundary faces
     115         [ +  + ]:     426426 :     for (std::size_t e=0; e<esuel.size()/4; ++e) {
     116                 :     426420 :       auto mark = e*4;
     117         [ +  + ]:    2132100 :       for (std::size_t f=0; f<4; ++f) {
     118         [ +  + ]:    1705680 :         if (esuel[mark+f] == -1) {
     119                 :            :           // extract triangle element connectivity from tetrahedron
     120         [ +  - ]:      46712 :           triinpoel.push_back( inpoel[ mark+tk::lpofa[f][0] ] );
     121         [ +  - ]:      46712 :           triinpoel.push_back( inpoel[ mark+tk::lpofa[f][1] ] );
     122         [ +  - ]:      46712 :           triinpoel.push_back( inpoel[ mark+tk::lpofa[f][2] ] );
     123                 :            :         }
     124                 :            :       }
     125                 :            :     }
     126                 :            : 
     127 [ +  - ][ +  - ]:          6 :     print.diagend( "done" );
     128 [ +  - ][ +  - ]:          6 :     times.emplace_back( "Generate surface mesh", t.dsec() );
     129                 :          6 :     t.zero();
     130                 :            :   }
     131                 :            : 
     132         [ +  + ]:         18 :   if (reorder) {
     133 [ +  - ][ +  - ]:          1 :     print.diagstart( "Reordering mesh nodes ..." );
     134                 :            : 
     135                 :            :     // If mesh has tetrahedra elements, reorder based on those
     136         [ +  - ]:          1 :     if (!mesh.tetinpoel().empty()) {
     137                 :            : 
     138                 :          1 :       auto& inpoel = mesh.tetinpoel();
     139 [ +  - ][ +  - ]:          2 :       const auto psup = tk::genPsup( inpoel, 4, tk::genEsup( inpoel, 4 ) );
     140         [ +  - ]:          2 :       auto map = tk::renumber( psup );
     141         [ +  - ]:          1 :       tk::remap( inpoel, map );
     142         [ +  - ]:          1 :       tk::remap( mesh.triinpoel(), map );
     143         [ +  - ]:          1 :       tk::remap( mesh.x(), map );
     144         [ +  - ]:          1 :       tk::remap( mesh.y(), map );
     145         [ +  - ]:          1 :       tk::remap( mesh.z(), map );
     146                 :            : 
     147                 :            :     // If mesh has no tetrahedra elements, reorder based on triangle mesh if any
     148         [ -  - ]:          0 :     } else if (!mesh.triinpoel().empty()) {
     149                 :            : 
     150                 :          0 :       auto& inpoel = mesh.triinpoel();
     151 [ -  - ][ -  - ]:          0 :       const auto psup = tk::genPsup( inpoel, 3, tk::genEsup( inpoel, 3 ) );
     152         [ -  - ]:          0 :       auto map = tk::renumber( psup );
     153         [ -  - ]:          0 :       tk::remap( inpoel, map );
     154         [ -  - ]:          0 :       tk::remap( mesh.x(), map );
     155         [ -  - ]:          0 :       tk::remap( mesh.y(), map );
     156         [ -  - ]:          0 :       tk::remap( mesh.z(), map );
     157                 :            :     }
     158                 :            : 
     159 [ +  - ][ +  - ]:          1 :     print.diagend( "done" );
     160 [ +  - ][ +  - ]:          1 :     times.emplace_back( "Reorder mesh", t.dsec() );
     161                 :          1 :     t.zero();
     162                 :            :   }
     163                 :            : 
     164 [ +  - ][ +  - ]:         18 :   print.diagstart( "Writing mesh to file '" + filename + "' ..." );
                 [ +  - ]
     165                 :            : 
     166         [ +  - ]:         18 :   const auto meshtype = pickOutput( filename );
     167                 :            : 
     168         [ +  + ]:         18 :   if (meshtype == MeshWriterType::GMSH)
     169 [ +  - ][ +  - ]:          4 :     GmshMeshWriter( filename ).writeMesh( mesh );
     170         [ +  + ]:         14 :   else if (meshtype == MeshWriterType::NETGEN)
     171 [ +  - ][ +  - ]:          4 :     NetgenMeshWriter( filename ).writeMesh( mesh );
     172         [ +  - ]:         10 :   else if (meshtype== MeshWriterType::EXODUSII)
     173 [ +  - ][ +  - ]:         10 :     ExodusIIMeshWriter( filename, ExoWriter::CREATE ).writeMesh( mesh );
     174                 :            : 
     175 [ +  - ][ +  - ]:         18 :   print.diagend( "done" );
     176 [ +  - ][ +  - ]:         18 :   times.emplace_back( "Write mesh to file", t.dsec() );
     177                 :            : 
     178                 :         36 :   return times;
     179                 :            : }
     180                 :            : 
     181                 :            : } // tk::