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