Chaste  Release::2018.1
GmshMeshReader.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF 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 
42 template<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();
72  ReadHeaders();
73 }
74 
75 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
77 {
78  CloseFiles();
79 }
80 
81 template<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 
94 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
96 {
97  mNodeFile.close();
98  mElementFile.close();
99  mFaceFile.close();
100 }
101 
102 template<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 
145 template<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 
163 template<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 
208 template<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 
250 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
252 {
253  return mNumElements;
254 }
255 
256 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
258 {
259  return mNumNodes;
260 }
261 
262 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
264 {
265  return mNumFaces;
266 }
267 
268 // LCOV_EXCL_START
269 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
271 {
273  //return mNumCableElements;
274 }
275 // LCOV_EXCL_STOP
276 
277 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
279 {
280  return mOrderOfElements;
281 }
282 
283 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
285 {
286  return mOrderOfBoundaryElements;
287 }
288 
289 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
291 {
292  return mNumElementAttributes;
293 }
294 
295 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
297 {
298  return mNumFaceAttributes;
299 }
300 
301 // LCOV_EXCL_START
302 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
304 {
306  //return mNumCableElementAttributes;
307 }
308 // LCOV_EXCL_STOP
309 
310 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
312 {
313  CloseFiles();
314 
315  OpenFiles();
316  ReadHeaders();
317 }
318 
319 template<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 
340 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
342 {
343  return std::vector<double>(0);
344 }
345 
346 template<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
401 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
403 {
405 }
406 // LCOV_EXCL_STOP
407 
408 template<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
463 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
464 std::vector<double> GmshMeshReader<ELEMENT_DIM, SPACE_DIM>::GetNode(unsigned index)
465 {
467 }
468 // LCOV_EXCL_STOP
469 
470 // LCOV_EXCL_START
471 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
473 {
475 }
476 // LCOV_EXCL_STOP
477 
478 // LCOV_EXCL_START
479 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
481 {
483 }
484 // LCOV_EXCL_STOP
485 
486 // Explicit instantiation
487 template class GmshMeshReader<0,1>;
488 template class GmshMeshReader<1,1>;
489 template class GmshMeshReader<1,2>;
490 template class GmshMeshReader<1,3>;
491 template class GmshMeshReader<2,2>;
492 template class GmshMeshReader<2,3>;
493 template class GmshMeshReader<3,3>;
unsigned GetNumElements() const
unsigned GetOrderOfBoundaryElements()
unsigned GetNumFaceAttributes() const
unsigned mOrderOfBoundaryElements
ElementData GetFaceData(unsigned index)
ElementData GetNextCableElementData()
unsigned mNodesPerBoundaryElement
unsigned GetOrderOfElements()
ElementData GetNextFaceData()
#define EXCEPTION(message)
Definition: Exception.hpp:143
#define NEVER_REACHED
Definition: Exception.hpp:206
ElementData GetNextElementData()
std::vector< unsigned > NodeIndices
unsigned GetNumFaces() const
std::vector< double > GetNodeAttributes()
ElementData GetElementData(unsigned index)
GmshMeshReader(std::string pathBaseName, unsigned orderOfElements=1, unsigned orderOfBoundaryElements=1)
std::vector< double > GetNode(unsigned index)
unsigned GetNumCableElements() const
unsigned GetNumElementAttributes() const
unsigned mOrderOfElements
unsigned GetNumCableElementAttributes() const
unsigned GetNumNodes() const
unsigned mNodesPerElement
std::vector< double > GetNextNode()