Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
GmshMeshReader.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35#include <cassert>
36#include <sstream>
37#include <iostream>
38
39#include "GmshMeshReader.hpp"
40#include "Exception.hpp"
41
42template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
44 unsigned orderOfElements,
45 unsigned orderOfBoundaryElements) :
46 mFileName(pathBaseName),
47 mOrderOfElements(orderOfElements),
48 mOrderOfBoundaryElements(orderOfBoundaryElements)
49{
50 // Only linear and quadratic elements
51 assert(mOrderOfElements==1 || mOrderOfElements==2);
52
53 if (mOrderOfElements==1)
54 {
55 mNodesPerElement = ELEMENT_DIM+1;
56 }
57 else
58 {
59 mNodesPerElement = (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2;
60 }
61
63 {
64 mNodesPerBoundaryElement = ELEMENT_DIM;
65 }
66 else
67 {
68 mNodesPerBoundaryElement = ELEMENT_DIM*(ELEMENT_DIM+1)/2;
69 }
70
71 OpenFiles();
73}
74
75template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
80
81template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
83{
84 // Open mesh file
85 mNodeFile.open(mFileName.c_str());
86 mElementFile.open(mFileName.c_str());
87 mFaceFile.open(mFileName.c_str());
88 if (!mNodeFile.is_open() || !mElementFile.is_open() || !mFaceFile.is_open() )
89 {
90 EXCEPTION("Could not open data file: " + mFileName);
91 }
92}
93
94template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
96{
97 mNodeFile.close();
98 mElementFile.close();
99 mFaceFile.close();
100}
101
102template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
104{
105 /*
106 * Read mesh format information from the file header
107 */
108 std::string this_line;
109 getline(mNodeFile, this_line);
110
111 assert(this_line == "$MeshFormat");
112
113 //Read the version no.
114 getline(mNodeFile, this_line);
115 std::stringstream line(this_line);
116
117 line >> mVersionNumber >> mFileType >> mDataSize;
118
119 if (mVersionNumber != 2.2)
120 {
121 EXCEPTION("Only .msh version 2.2 files are supported.");
122 }
123 assert(mFileType == 0);
124
125 //Check mesh format close string
126 getline(mNodeFile, this_line);
127 assert(this_line == "$EndMeshFormat");
128
129 ReadNodeHeader();
130 ReadElementHeader(); // This reads the total number of elements in the file into mTotalNumElementsAndFaces
131 ReadFaceHeader();
132
133 if (mTotalNumElementsAndFaces != mNumElements + mNumFaces)
134 {
135 EXCEPTION("Unrecognised element types present in the .msh file: check mesh generation settings in gmsh.");
136 // If you hit this, then the .msh file lists elements that are not simply
137 // in 2D:
138 // linear/quadratic lines for faces, triangles for elements
139 // in 3D:
140 // linear/quadratic traingles for faces, tets for elements.
141 // which probably means something funny happened in the mesh generation.
142 }
143}
144
145template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
147{
148 // Search for the start of the node section
149 std::string this_line;
150 std::stringstream line(this_line);
151
152 while (this_line != "$Nodes")
153 {
154 getline(mNodeFile, this_line);
155 }
156 getline(mNodeFile, this_line);
157
158 line.clear();
159 line.str(this_line);
160 line >> mNumNodes; // mNodesFile should now be pointing at the start of the node lines in the file.
161}
162
163template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
165{
166 // Search for the start of the elements section
167 std::string this_line;
168 std::stringstream line(this_line);
169
170 getline(mElementFile, this_line);
171 while (this_line != "$Elements")
172 {
173 getline(mElementFile, this_line);
174 }
175 getline(mElementFile, this_line); //Throw away the line that says the number of elements specified in the file
176 line.clear();
177 line.str(this_line);
178 line >> mTotalNumElementsAndFaces;
179 int ele_start = mElementFile.tellg(); //Pointer to the start of the element block.
180
181 mNumElements = 0u;
182 getline(mElementFile, this_line);
183
184 while (this_line != "$EndElements")
185 {
186 line.clear();
187 line.str(this_line);
188
189 unsigned ele_index;
190 unsigned ele_type;
191 line >> ele_index >> ele_type;
192
193 if (ELEMENT_DIM == 2 && (ele_type == GmshTypes::TRIANGLE || ele_type == GmshTypes::QUADRATIC_TRIANGLE))
194 {
195 mNumElements++;
196 }
197 else if (ELEMENT_DIM == 3 && (ele_type == GmshTypes::TETRAHEDRON || ele_type == GmshTypes::QUADRATIC_TETRAHEDRON))
198 {
199 mNumElements++;
200 }
201
202 getline(mElementFile, this_line);
203 }
204
205 mElementFile.seekg(ele_start); //mElementFile should now be pointing at the start of the node lines in the file.
206}
207
208template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
210{
211 // Search for the start of the elements section
212 std::string this_line;
213 std::stringstream line(this_line);
214
215 getline(mFaceFile, this_line);
216 while (this_line != "$Elements")
217 {
218 getline(mFaceFile, this_line);
219 }
220 getline(mFaceFile, this_line); // Throw away the line that says the number of elements specified in the file
221 int face_start = mFaceFile.tellg(); // Pointer to the start of the element block.
222
223 mNumFaces = 0u;
224 getline(mFaceFile, this_line);
225
226 while (this_line != "$EndElements")
227 {
228 line.clear();
229 line.str(this_line);
230
231 unsigned ele_index;
232 unsigned ele_type;
233 line >> ele_index >> ele_type;
234
235 if (ELEMENT_DIM == 2 && (ele_type == GmshTypes::LINE || ele_type == GmshTypes::QUADRATIC_LINE))
236 {
237 mNumFaces++;
238 }
239 else if (ELEMENT_DIM == 3 && (ele_type == GmshTypes::TRIANGLE || ele_type == GmshTypes::QUADRATIC_TRIANGLE))
240 {
241 mNumFaces++;
242 }
243
244 getline(mFaceFile, this_line);
245 }
246
247 mFaceFile.seekg(face_start); //mFacesFile should now be pointing at the start of the node lines in the file.
248}
249
250template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
252{
253 return mNumElements;
254}
255
256template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
258{
259 return mNumNodes;
260}
261
262template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
264{
265 return mNumFaces;
266}
267
268// LCOV_EXCL_START
269template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
271{
273 //return mNumCableElements;
274}
275// LCOV_EXCL_STOP
276
277template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
279{
280 return mOrderOfElements;
281}
282
283template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
285{
286 return mOrderOfBoundaryElements;
287}
288
289template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
291{
292 return mNumElementAttributes;
293}
294
295template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
297{
298 return mNumFaceAttributes;
299}
300
301// LCOV_EXCL_START
302template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
304{
306 //return mNumCableElementAttributes;
307}
308// LCOV_EXCL_STOP
309
310template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
312{
313 CloseFiles();
314
315 OpenFiles();
316 ReadHeaders();
317}
318
319template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
321{
322 std::vector<double> ret_coords(SPACE_DIM);
323
324 std::string this_line;
325 getline(mNodeFile, this_line);
326
327 std::stringstream line(this_line);
328
329 unsigned node_index;
330 line >> node_index >> ret_coords[0] >> ret_coords[1];
331
332 if (SPACE_DIM == 3)
333 {
334 line >> ret_coords[2];
335 }
336
337 return ret_coords;
338}
339
340template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
342{
343 return std::vector<double>(0);
344}
345
346template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
348{
349 ElementData element_data;
350 element_data.NodeIndices.resize(mNodesPerElement);
351 element_data.AttributeValue = 0.0; // If an attribute is not read this stays as zero, otherwise overwritten.
352
353 std::string this_line;
354 std::stringstream line(this_line);
355
356 unsigned ele_index;
357 unsigned ele_type;
358 unsigned ele_attributes = 0;
359 bool volume_element_found = false;
360
361 while (!volume_element_found)
362 {
363 getline(mElementFile, this_line);
364 line.clear();
365 line.str(this_line);
366
367 line >> ele_index >> ele_type >> ele_attributes;
368
369 if ((ELEMENT_DIM == 2 && (ele_type == GmshTypes::TRIANGLE || ele_type == GmshTypes::QUADRATIC_TRIANGLE) ) ||
370 (ELEMENT_DIM == 3 && (ele_type == GmshTypes::TETRAHEDRON || ele_type == GmshTypes::QUADRATIC_TETRAHEDRON)))
371 {
372 volume_element_found = true;
373 }
374 }
375
376 //Gmsh can have arbitrary numbers of element attributes, but Chaste only handles one. We pick the first attribute and throw
377 //away the remainder.
378 if (ele_attributes > 0)
379 {
380 mNumElementAttributes = 1u;
381 line >> element_data.AttributeValue;
382 }
383 unsigned unused_attr;
384 for (unsigned attr_index = 0; attr_index < (ele_attributes-1); ++attr_index)
385 {
386 line >> unused_attr;
387 }
388
389 //Read the node indices
390 for (unsigned node_index = 0; node_index < mNodesPerElement; ++node_index)
391 {
392 line >> element_data.NodeIndices[node_index];
393
394 element_data.NodeIndices[node_index]--; //Gmsh *always* indexes from 1, we index from 0
395 }
396
397 return element_data;
398}
399
400// LCOV_EXCL_START
401template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
406// LCOV_EXCL_STOP
407
408template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
410{
411 ElementData face_data;
412 face_data.NodeIndices.resize(mNodesPerBoundaryElement);
413 face_data.AttributeValue = 0.0; // If an attribute is not read this stays as zero, otherwise overwritten.
414
415 std::string this_line;
416 std::stringstream line(this_line);
417
418 unsigned face_index;
419 unsigned face_type;
420 unsigned face_attributes=0;
421 bool surface_element_found = false;
422
423 while (!surface_element_found)
424 {
425 getline(mFaceFile, this_line);
426 line.clear();
427 line.str(this_line);
428
429 line >> face_index >> face_type >> face_attributes;
430
431 if ((ELEMENT_DIM == 2 && (face_type == GmshTypes::LINE || face_type == GmshTypes::QUADRATIC_LINE) ) ||
432 (ELEMENT_DIM == 3 && (face_type == GmshTypes::TRIANGLE || face_type == GmshTypes::QUADRATIC_TRIANGLE)))
433 {
434 surface_element_found = true;
435 }
436 }
437
438 //Gmsh can have arbitrary numbers of element attributes, but Chaste only handles one. We pick the first attribute and throw
439 //away the remainder.
440 if (face_attributes > 0)
441 {
442 mNumFaceAttributes = 1u;
443 line >> face_data.AttributeValue;
444 }
445 unsigned unused_attr;
446 for (unsigned attr_index = 0; attr_index < (face_attributes-1); ++attr_index)
447 {
448 line >> unused_attr;
449 }
450
451 //Read the node indices
452 for (unsigned node_index = 0; node_index < mNodesPerBoundaryElement; ++node_index)
453 {
454 line >> face_data.NodeIndices[node_index];
455
456 face_data.NodeIndices[node_index]--; //Gmsh *always* indexes from 1, we index from 0
457 }
458
459 return face_data;
460}
461
462// LCOV_EXCL_START
463template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
464std::vector<double> GmshMeshReader<ELEMENT_DIM, SPACE_DIM>::GetNode(unsigned index)
465{
467}
468// LCOV_EXCL_STOP
469
470// LCOV_EXCL_START
471template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
476// LCOV_EXCL_STOP
477
478// LCOV_EXCL_START
479template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
484// LCOV_EXCL_STOP
485
486// Explicit instantiation
487template class GmshMeshReader<0,1>;
488template class GmshMeshReader<1,1>;
489template class GmshMeshReader<1,2>;
490template class GmshMeshReader<1,3>;
491template class GmshMeshReader<2,2>;
492template class GmshMeshReader<2,3>;
493template class GmshMeshReader<3,3>;
#define EXCEPTION(message)
#define NEVER_REACHED
ElementData GetNextFaceData()
unsigned mOrderOfElements
unsigned GetNumElements() const
GmshMeshReader(std::string pathBaseName, unsigned orderOfElements=1, unsigned orderOfBoundaryElements=1)
unsigned GetNumCableElements() const
ElementData GetNextCableElementData()
unsigned GetOrderOfElements()
unsigned mOrderOfBoundaryElements
std::vector< double > GetNode(unsigned index)
unsigned mNodesPerElement
unsigned GetNumNodes() const
ElementData GetNextElementData()
unsigned mNodesPerBoundaryElement
unsigned GetNumElementAttributes() const
ElementData GetElementData(unsigned index)
unsigned GetNumCableElementAttributes() const
unsigned GetNumFaceAttributes() const
unsigned GetNumFaces() const
std::vector< double > GetNodeAttributes()
ElementData GetFaceData(unsigned index)
unsigned GetOrderOfBoundaryElements()
std::vector< double > GetNextNode()
std::vector< unsigned > NodeIndices