Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/IO/GmshMeshReader.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 Gmsh mesh reader class definition
10 : : \details Gmsh mesh reader class definition. Currently, this class supports
11 : : line, triangle, tetrahedron, and point Gmsh element types.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <limits>
16 : : #include <array>
17 : : #include <cmath>
18 : : #include <cstddef>
19 : : #include <istream>
20 : : #include <string>
21 : : #include <utility>
22 : : #include <vector>
23 : :
24 : : #include "UnsMesh.hpp"
25 : : #include "GmshMeshReader.hpp"
26 : : #include "GmshMeshIO.hpp"
27 : : #include "Reorder.hpp"
28 : : #include "PrintUtil.hpp"
29 : :
30 : : using tk::GmshMeshReader;
31 : :
32 : : void
33 : 7 : GmshMeshReader::readMesh( UnsMesh& mesh )
34 : : // *****************************************************************************
35 : : // Public interface for read a Gmsh mesh from file
36 : : //! \param[in] mesh Unstructured mesh object
37 : : // *****************************************************************************
38 : : {
39 : : // Read in mandatory "$MeshFormat" section
40 : 7 : readMeshFormat();
41 : :
42 : : // Keep reading in sections until end of file. These sections can be in
43 : : // arbitrary order, hence a while loop.
44 [ + + ]: 28 : while ( !m_inFile.eof() ) {
45 : 21 : std::string s;
46 [ + - ]: 21 : getline( m_inFile, s );
47 [ + + ]: 21 : if ( s == "$Nodes" )
48 [ + - ]: 7 : readNodes( mesh );
49 [ + + ]: 14 : else if ( s == "$Elements" )
50 [ + - ]: 7 : readElements( mesh );
51 : 21 : }
52 : 7 : }
53 : :
54 : : void
55 : 7 : GmshMeshReader::readMeshFormat()
56 : : // *****************************************************************************
57 : : // Read mandatory "$MeshFormat--$EndMeshFormat" section
58 : : // *****************************************************************************
59 : : {
60 : : using tk::operator<<;
61 : 7 : std::string s;
62 : :
63 : : // Read in beginning of header: $MeshFormat
64 [ + - ]: 7 : getline( m_inFile, s );
65 [ - + ][ - - ]: 7 : ErrChk( s == "$MeshFormat",
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
66 : : std::string("Unsupported mesh format '") + s + "' in file " +
67 : : m_filename );
68 : :
69 : : // Read in "version-number file-type data-size"
70 : : int type;
71 [ + - ][ + - ]: 7 : m_inFile >> m_version >> type >> m_datasize;
[ + - ]
72 [ + + ]: 7 : if (type == 0 )
73 : 4 : m_type = GmshFileType::ASCII;
74 [ + - ]: 3 : else if (type == 1 )
75 : 3 : m_type = GmshFileType::BINARY;
76 : :
77 [ - + ][ - - ]: 7 : ErrChk( ( fabs(m_version-2.2) < std::numeric_limits<tk::real>::epsilon() ||
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
78 : : fabs(m_version-2.0) < std::numeric_limits<tk::real>::epsilon() ),
79 : : std::string("Unsupported mesh version '") << m_version <<
80 : : "' in file " << m_filename );
81 : :
82 [ + + ][ - + ]: 7 : ErrChk( ( m_type == GmshFileType::ASCII || m_type == GmshFileType::BINARY ),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
83 : : std::string("Unsupported mesh type '") << m_type << "' in file " <<
84 : : m_filename );
85 : :
86 [ - + ][ - - ]: 7 : ErrChk( m_datasize == sizeof(tk::real),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
87 : : std::string("Unsupported mesh datasize '") << m_datasize <<
88 : : "' in file " << m_filename );
89 : :
90 [ + - ]: 7 : getline( m_inFile, s ); // finish reading the line
91 : :
92 : : // if file is binary, binary-read in binary "one"
93 [ + - ][ + + ]: 7 : if ( isBinary() ) {
94 : : int one;
95 [ + - ]: 3 : m_inFile.read( reinterpret_cast<char*>(&one), sizeof(int) );
96 : : #ifdef __bg__
97 : : one = tk::swap_endian< int >( one );
98 : : #endif
99 [ - + ][ - - ]: 3 : ErrChk( one == 1, "Endianness does not match in file " + m_filename );
[ - - ][ - - ]
100 [ + - ]: 3 : getline( m_inFile, s ); // finish reading the line
101 : : }
102 : :
103 : : // Read in end of header: $EndMeshFormat
104 [ + - ]: 7 : getline( m_inFile, s );
105 [ - + ][ - - ]: 7 : ErrChk( s == "$EndMeshFormat",
[ - - ][ - - ]
106 : : "'$EndMeshFormat' keyword is missing in file " + m_filename );
107 : 7 : }
108 : :
109 : : void
110 : 7 : GmshMeshReader::readNodes( UnsMesh& mesh )
111 : : // *****************************************************************************
112 : : // Read "$Nodes--$EndNodes" section
113 : : //! \param[in] mesh Unstructured mesh object
114 : : // *****************************************************************************
115 : : {
116 : : // Read in number of nodes in this node set
117 : : std::size_t nnode;
118 [ + - ]: 7 : m_inFile >> nnode;
119 [ - + ][ - - ]: 7 : ErrChk( nnode > 0,
[ - - ][ - - ]
120 : : "Number of nodes must be greater than zero in file " + m_filename );
121 : 7 : std::string s;
122 [ + - ][ + + ]: 7 : if (isBinary()) getline( m_inFile, s ); // finish reading the line
[ + - ]
123 : :
124 : : // Read in node ids and coordinates: node-number x-coord y-coord z-coord
125 [ + + ]: 105 : for ( std::size_t i=0; i<nnode; ++i ) {
126 : : int id;
127 : : std::array< tk::real, 3 > coord;
128 : :
129 [ + - ][ + + ]: 98 : if (isASCII()) {
130 [ + - ][ + - ]: 56 : m_inFile >> id >> coord[0] >> coord[1] >> coord[2];
[ + - ][ + - ]
131 : : } else {
132 [ + - ]: 42 : m_inFile.read( reinterpret_cast<char*>(&id), sizeof(int) );
133 : : #ifdef __bg__
134 : : id = tk::swap_endian< int >( id );
135 : : #endif
136 [ + - ]: 42 : m_inFile.read( reinterpret_cast<char*>(coord.data()), 3*sizeof(double) );
137 : : #ifdef __bg__
138 : : coord[0] = tk::swap_endian< double >( coord[0] );
139 : : coord[1] = tk::swap_endian< double >( coord[1] );
140 : : coord[2] = tk::swap_endian< double >( coord[2] );
141 : : #endif
142 : : }
143 : :
144 [ + - ]: 98 : mesh.x().push_back( coord[0] );
145 [ + - ]: 98 : mesh.y().push_back( coord[1] );
146 [ + - ]: 98 : mesh.z().push_back( coord[2] );
147 : : }
148 [ + - ]: 7 : getline( m_inFile, s ); // finish reading the last line
149 : :
150 : : // Read in end of header: $EndNodes
151 [ + - ]: 7 : getline( m_inFile, s );
152 [ - + ][ - - ]: 7 : ErrChk( s == "$EndNodes",
[ - - ][ - - ]
153 : : "'$EndNodes' keyword is missing in file" + m_filename );
154 : 7 : }
155 : :
156 : : void
157 : 7 : GmshMeshReader::readElements( UnsMesh& mesh )
158 : : // *****************************************************************************
159 : : // Read "$Elements--$EndElements" section
160 : : //! \param[in] mesh Unstructured mesh object
161 : : // *****************************************************************************
162 : : {
163 : : using tk::operator<<;
164 : :
165 : 7 : std::string s;
166 : :
167 : : // Read in number of elements in this element set
168 : : int nel;
169 [ + - ]: 7 : m_inFile >> nel;
170 [ - + ][ - - ]: 7 : ErrChk( nel > 0, "Number of elements must be greater than zero in file " +
[ - - ][ - - ]
171 : : m_filename );
172 [ + - ]: 7 : getline( m_inFile, s ); // finish reading the last line
173 : :
174 : : // Read in element ids, tags, and element connectivity (node list)
175 : 7 : int n=1;
176 [ + + ]: 180 : for (int i=0; i<nel; i+=n) {
177 : : int id, elmtype, ntags;
178 : :
179 [ + - ][ + + ]: 173 : if (isASCII()) {
180 : : // elm-number elm-type number-of-tags < tag > ... node-number-list
181 [ + - ][ + - ]: 168 : m_inFile >> id >> elmtype >> ntags;
[ + - ]
182 : : } else {
183 : : // elm-type num-of-elm-follow number-of-tags
184 [ + - ]: 5 : m_inFile.read( reinterpret_cast<char*>(&elmtype), sizeof(int) );
185 [ + - ]: 5 : m_inFile.read( reinterpret_cast<char*>(&n), sizeof(int) );
186 [ + - ]: 5 : m_inFile.read( reinterpret_cast<char*>(&ntags), sizeof(int) );
187 : : #ifdef __bg__
188 : : elmtype = tk::swap_endian< int >( elmtype );
189 : : n = tk::swap_endian< int >( n );
190 : : ntags = tk::swap_endian< int >( ntags );
191 : : #endif
192 : : }
193 : :
194 : : // Find element type, throw exception if not supported
195 [ + - ]: 173 : const auto it = m_elemNodes.find( elmtype );
196 [ - + ][ - - ]: 173 : ErrChk( it != m_elemNodes.end(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
197 : : std::string("Unsupported element type ") << elmtype <<
198 : : " in mesh file: " << m_filename );
199 : :
200 [ + + ]: 461 : for (int e=0; e<n; ++e) {
201 : : // Read element id if binary
202 [ + - ][ + + ]: 288 : if (isBinary()) {
203 [ + - ]: 120 : m_inFile.read( reinterpret_cast<char*>(&id), sizeof(int) );
204 : : #ifdef __bg__
205 : : id = tk::swap_endian< int >( id );
206 : : #endif
207 : : }
208 : :
209 : : // Read and ignore element tags
210 [ + - ]: 288 : std::vector< int > tags( static_cast<std::size_t>(ntags), 0 );
211 [ + - ][ + + ]: 288 : if (isASCII()) {
212 [ + + ]: 480 : for (std::size_t j=0; j<static_cast<std::size_t>(ntags); j++)
213 [ + - ]: 312 : m_inFile >> tags[j];
214 : : } else {
215 [ + - ]: 120 : m_inFile.read(
216 : 120 : reinterpret_cast<char*>(tags.data()),
217 : : static_cast<std::streamsize>(
218 : 120 : static_cast<std::size_t>(ntags) * sizeof(int) ) );
219 : : #ifdef __bg__
220 : : for (auto& t : tags) t = tk::swap_endian< int >( t );
221 : : #endif
222 : : }
223 : :
224 : : // Read and add element node list (i.e. connectivity)
225 : 288 : std::size_t nnode = static_cast< std::size_t >( it->second );
226 [ + - ]: 288 : std::vector< std::size_t > nodes( nnode, 0 );
227 [ + - ][ + + ]: 288 : if (isASCII()) {
228 [ + + ]: 768 : for (std::size_t j=0; j<nnode; j++)
229 [ + - ]: 600 : m_inFile >> nodes[j];
230 : : } else {
231 [ + - ]: 120 : std::vector< int > nds( nnode, 0 );
232 [ + - ]: 120 : m_inFile.read(
233 : 120 : reinterpret_cast< char* >( nds.data() ),
234 : 120 : static_cast< std::streamsize >( nnode * sizeof(int) ) );
235 : : #ifdef __bg__
236 : : for (auto& j : nds) j = tk::swap_endian< int >( j );
237 : : #endif
238 [ + + ]: 552 : for (std::size_t j=0; j<nnode; j++)
239 : 432 : nodes[j] = static_cast< std::size_t >( nds[j] );
240 : 120 : }
241 : : // Put in element connectivity for different types of elements
242 [ + + ][ - - ]: 288 : switch ( elmtype ) {
243 : 120 : case GmshElemType::TRI:
244 [ + - ][ + + ]: 480 : for (const auto& j : nodes) mesh.triinpoel().push_back( j );
245 : 120 : break;
246 : 168 : case GmshElemType::TET:
247 [ + - ][ + + ]: 840 : for (const auto& j : nodes) mesh.tetinpoel().push_back( j );
248 : 168 : break;
249 : 0 : case GmshElemType::PNT:
250 : 0 : break; // ignore 1-node 'point element' type
251 [ - - ][ - - ]: 0 : default: Throw( std::string("Unsupported element type ") << elmtype <<
[ - - ][ - - ]
[ - - ][ - - ]
252 : : " in mesh file: " << m_filename );
253 : : }
254 : 288 : }
255 : : }
256 [ + - ]: 7 : getline( m_inFile, s ); // finish reading the last line
257 : :
258 : : // Shift node IDs to start from zero (gmsh likes one-based node ids)
259 [ + - ]: 7 : shiftToZero( mesh.triinpoel() );
260 [ + - ]: 7 : shiftToZero( mesh.tetinpoel() );
261 : :
262 : : // Read in end of header: $EndNodes
263 [ + - ]: 7 : getline( m_inFile, s );
264 [ - + ][ - - ]: 7 : ErrChk( s == "$EndElements",
[ - - ][ - - ]
265 : : "'$EndElements' keyword is missing in file" + m_filename );
266 : 7 : }
|
