Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
InputGmsh.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: InputGmsh.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: GMSH converter.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "MeshElements.h"
41 #include "InputGmsh.h"
42 
43 namespace Nektar
44 {
45  namespace Utilities
46  {
47  ModuleKey InputGmsh::className =
48  GetModuleFactory().RegisterCreatorFunction(
49  ModuleKey(eInputModule, "msh"), InputGmsh::create,
50  "Reads Gmsh msh file.");
51 
52  std::map<unsigned int, ElmtConfig> InputGmsh::elmMap =
53  InputGmsh::GenElmMap();
54 
55 
56  /**
57  * @brief Set up InputGmsh object.
58  *
59  */
60  InputGmsh::InputGmsh(MeshSharedPtr m) : InputModule(m)
61  {
62 
63  }
64 
65  InputGmsh::~InputGmsh()
66  {
67 
68  }
69 
70  /**
71  * Gmsh file contains a list of nodes and their coordinates, along with
72  * a list of elements and those nodes which define them. We read in and
73  * store the list of nodes in #m_node and store the list of elements in
74  * #m_element. Each new element is supplied with a list of entries from
75  * #m_node which defines the element. Finally some mesh statistics are
76  * printed.
77  *
78  * @param pFilename Filename of Gmsh file to read.
79  */
80  void InputGmsh::Process()
81  {
82  // Open the file stream.
83  OpenStream();
84 
85  m_mesh->m_expDim = 0;
86  m_mesh->m_spaceDim = 0;
87  string line;
88  int nVertices = 0;
89  int nEntities = 0;
90  int elm_type = 0;
91  int prevId = -1;
92  map<unsigned int, ElmtConfig>::iterator it;
93 
94  if (m_mesh->m_verbose)
95  {
96  cout << "InputGmsh: Start reading file..." << endl;
97  }
98 
99  while (!m_mshFile.eof())
100  {
101  getline(m_mshFile, line);
102  stringstream s(line);
103  string word;
104  s >> word;
105 
106  // Process nodes.
107  if (word == "$Nodes")
108  {
109  getline(m_mshFile, line);
110  stringstream s(line);
111  s >> nVertices;
112  int id = 0;
113  for (int i = 0; i < nVertices; ++i)
114  {
115  getline(m_mshFile, line);
116  stringstream st(line);
117  double x = 0, y = 0, z = 0;
118  st >> id >> x >> y >> z;
119 
120  if ((x * x) > 0.000001 && m_mesh->m_spaceDim < 1)
121  {
122  m_mesh->m_spaceDim = 1;
123  }
124  if ((y * y) > 0.000001 && m_mesh->m_spaceDim < 2)
125  {
126  m_mesh->m_spaceDim = 2;
127  }
128  if ((z * z) > 0.000001 && m_mesh->m_spaceDim < 3)
129  {
130  m_mesh->m_spaceDim = 3;
131  }
132 
133  id -= 1; // counter starts at 0
134 
135  if (id-prevId != 1)
136  {
137  cerr << "Gmsh vertex ids should be contiguous" << endl;
138  abort();
139  }
140  prevId = id;
141  m_mesh->m_node.push_back(boost::shared_ptr<Node>(new Node(id, x, y, z)));
142  }
143  }
144  // Process elements
145  else if (word == "$Elements")
146  {
147  getline(m_mshFile, line);
148  stringstream s(line);
149  s >> nEntities;
150  for (int i = 0; i < nEntities; ++i)
151  {
152  getline(m_mshFile, line);
153  stringstream st(line);
154  int id = 0, num_tag = 0, num_nodes = 0;
155 
156  st >> id >> elm_type >> num_tag;
157  id -= 1; // counter starts at 0
158 
159  it = elmMap.find(elm_type);
160  if (it == elmMap.end())
161  {
162  cerr << "Error: element type " << elm_type
163  << " not supported" << endl;
164  abort();
165  }
166 
167  // Read element tags
168  vector<int> tags;
169  for (int j = 0; j < num_tag; ++j)
170  {
171  int tag = 0;
172  st >> tag;
173  tags.push_back(tag);
174  }
175  tags.resize(1);
176 
177  // Read element node list
178  vector<NodeSharedPtr> nodeList;
179  num_nodes = GetNnodes(elm_type);
180  for (int k = 0; k < num_nodes; ++k)
181  {
182  int node = 0;
183  st >> node;
184  node -= 1; // counter starts at 0
185  nodeList.push_back(m_mesh->m_node[node]);
186  }
187 
188  // Prism nodes need re-ordering for Nektar++.
189  if (it->second.m_e == LibUtilities::ePrism)
190  {
191  // Mirror first in uv plane to swap around
192  // triangular faces
193  swap(nodeList[0], nodeList[3]);
194  swap(nodeList[1], nodeList[4]);
195  swap(nodeList[2], nodeList[5]);
196  // Reorder base points so that face/vertices map
197  // correctly.
198  swap(nodeList[4], nodeList[2]);
199 
200  if (it->second.m_order == 2)
201  {
202  vector<NodeSharedPtr> nodemap(18);
203 
204  // Vertices remain unchanged.
205  nodemap[ 0] = nodeList[ 0];
206  nodemap[ 1] = nodeList[ 1];
207  nodemap[ 2] = nodeList[ 2];
208  nodemap[ 3] = nodeList[ 3];
209  nodemap[ 4] = nodeList[ 4];
210  nodemap[ 5] = nodeList[ 5];
211  // Reorder edge nodes: first mirror in uv
212  // plane and then place in Nektar++ ordering.
213  nodemap[ 6] = nodeList[12];
214  nodemap[ 7] = nodeList[10];
215  nodemap[ 8] = nodeList[ 6];
216  nodemap[ 9] = nodeList[ 8];
217  nodemap[10] = nodeList[13];
218  nodemap[11] = nodeList[14];
219  nodemap[12] = nodeList[ 9];
220  nodemap[13] = nodeList[ 7];
221  nodemap[14] = nodeList[11];
222  // Face vertices remain unchanged.
223  nodemap[15] = nodeList[15];
224  nodemap[16] = nodeList[16];
225  nodemap[17] = nodeList[17];
226 
227  nodeList = nodemap;
228  }
229  else if (it->second.m_order > 2)
230  {
231  cerr << "Error: gmsh prisms only supported up "
232  << "to second order." << endl;
233  abort();
234  }
235  }
236 
237  // Create element
238  ElementSharedPtr E = GetElementFactory().
239  CreateInstance(it->second.m_e,it->second,nodeList,tags);
240 
241  // Determine mesh expansion dimension
242  if (E->GetDim() > m_mesh->m_expDim) {
243  m_mesh->m_expDim = E->GetDim();
244  }
245  m_mesh->m_element[E->GetDim()].push_back(E);
246  }
247  }
248  }
249  m_mshFile.close();
250 
251  // Process rest of mesh.
252  ProcessVertices ();
253  ProcessEdges ();
254  ProcessFaces ();
255  ProcessElements ();
256  ProcessComposites();
257  }
258 
259  /**
260  * For a given msh ID, return the corresponding number of nodes.
261  */
262  int InputGmsh::GetNnodes(unsigned int InputGmshEntity)
263  {
264  int nNodes;
265  map<unsigned int, ElmtConfig>::iterator it;
266 
267  it = elmMap.find(InputGmshEntity);
268 
269  if (it == elmMap.end())
270  {
271  cerr << "Unknown element type " << InputGmshEntity << endl;
272  abort();
273  }
274 
275  switch(it->second.m_e)
276  {
277  case LibUtilities::ePoint:
278  nNodes = Point:: GetNumNodes(it->second);
279  break;
280  case LibUtilities::eSegment:
281  nNodes = Line:: GetNumNodes(it->second);
282  break;
283  case LibUtilities::eTriangle:
284  nNodes = Triangle:: GetNumNodes(it->second);
285  break;
286  case LibUtilities::eQuadrilateral:
287  nNodes = Quadrilateral::GetNumNodes(it->second);
288  break;
289  case LibUtilities::eTetrahedron:
290  nNodes = Tetrahedron:: GetNumNodes(it->second);
291  break;
292  case LibUtilities::ePyramid:
293  nNodes = Pyramid:: GetNumNodes(it->second);
294  break;
295  case LibUtilities::ePrism:
296  nNodes = Prism:: GetNumNodes(it->second);
297  break;
298  case LibUtilities::eHexahedron:
299  nNodes = Hexahedron:: GetNumNodes(it->second);
300  break;
301  default:
302  cerr << "Unknown element type!" << endl;
303  abort();
304  break;
305  }
306 
307  return nNodes;
308  }
309 
310  /*
311  * @brief Populate the element map #elmMap.
312  *
313  * This function primarily populates the element mapping #elmMap,
314  * which takes a msh ID used by Gmsh and translates to element type,
315  * element order and whether the element is incomplete (i.e. whether
316  * it contains solely boundary nodes, or just face nodes). Note that
317  * some of these elements, such as prisms of order >= 3, cannot yet be
318  * generated by Gmsh.
319  */
320  std::map<unsigned int, ElmtConfig> InputGmsh::GenElmMap()
321  {
322  using namespace LibUtilities;
323  std::map<unsigned int, ElmtConfig> tmp;
324 
325  // Elmt type, order, face, volume
326  tmp[ 1] = ElmtConfig(eSegment, 1, true, true);
327  tmp[ 2] = ElmtConfig(eTriangle, 1, true, true);
328  tmp[ 3] = ElmtConfig(eQuadrilateral, 1, true, true);
329  tmp[ 4] = ElmtConfig(eTetrahedron, 1, true, true);
330  tmp[ 5] = ElmtConfig(eHexahedron, 1, true, true);
331  tmp[ 6] = ElmtConfig(ePrism, 1, true, true);
332  tmp[ 7] = ElmtConfig(ePyramid, 1, true, true);
333  tmp[ 8] = ElmtConfig(eSegment, 2, true, true);
334  tmp[ 9] = ElmtConfig(eTriangle, 2, true, true);
335  tmp[ 10] = ElmtConfig(eQuadrilateral, 2, true, true);
336  tmp[ 11] = ElmtConfig(eTetrahedron, 2, true, true);
337  tmp[ 12] = ElmtConfig(eHexahedron, 2, true, true);
338  tmp[ 13] = ElmtConfig(ePrism, 2, true, true);
339  tmp[ 15] = ElmtConfig(ePoint, 1, true, false);
340  tmp[ 16] = ElmtConfig(eQuadrilateral, 2, false, false);
341  tmp[ 17] = ElmtConfig(eHexahedron, 2, false, false);
342  tmp[ 18] = ElmtConfig(ePrism, 2, false, false);
343  tmp[ 20] = ElmtConfig(eTriangle, 3, false, false);
344  tmp[ 21] = ElmtConfig(eTriangle, 3, true, false);
345  tmp[ 22] = ElmtConfig(eTriangle, 4, false, false);
346  tmp[ 23] = ElmtConfig(eTriangle, 4, true, false);
347  tmp[ 24] = ElmtConfig(eTriangle, 5, false, false);
348  tmp[ 25] = ElmtConfig(eTriangle, 5, true, false);
349  tmp[ 26] = ElmtConfig(eSegment, 3, true, false);
350  tmp[ 27] = ElmtConfig(eSegment, 4, true, false);
351  tmp[ 28] = ElmtConfig(eSegment, 5, true, false);
352  tmp[ 29] = ElmtConfig(eTetrahedron, 3, true, true);
353  tmp[ 30] = ElmtConfig(eTetrahedron, 4, true, true);
354  tmp[ 31] = ElmtConfig(eTetrahedron, 5, true, true);
355  tmp[ 32] = ElmtConfig(eTetrahedron, 4, true, false);
356  tmp[ 33] = ElmtConfig(eTetrahedron, 5, true, false);
357  tmp[ 36] = ElmtConfig(eQuadrilateral, 3, true, false);
358  tmp[ 37] = ElmtConfig(eQuadrilateral, 4, true, false);
359  tmp[ 38] = ElmtConfig(eQuadrilateral, 5, true, false);
360  tmp[ 39] = ElmtConfig(eQuadrilateral, 3, false, false);
361  tmp[ 40] = ElmtConfig(eQuadrilateral, 4, false, false);
362  tmp[ 41] = ElmtConfig(eQuadrilateral, 5, false, false);
363  tmp[ 42] = ElmtConfig(eTriangle, 6, true, false);
364  tmp[ 43] = ElmtConfig(eTriangle, 7, true, false);
365  tmp[ 44] = ElmtConfig(eTriangle, 8, true, false);
366  tmp[ 45] = ElmtConfig(eTriangle, 9, true, false);
367  tmp[ 46] = ElmtConfig(eTriangle, 10, true, false);
368  tmp[ 47] = ElmtConfig(eQuadrilateral, 6, true, false);
369  tmp[ 48] = ElmtConfig(eQuadrilateral, 7, true, false);
370  tmp[ 49] = ElmtConfig(eQuadrilateral, 8, true, false);
371  tmp[ 50] = ElmtConfig(eQuadrilateral, 9, true, false);
372  tmp[ 51] = ElmtConfig(eQuadrilateral, 10, true, false);
373  tmp[ 52] = ElmtConfig(eTriangle, 6, false, false);
374  tmp[ 53] = ElmtConfig(eTriangle, 7, false, false);
375  tmp[ 54] = ElmtConfig(eTriangle, 8, false, false);
376  tmp[ 55] = ElmtConfig(eTriangle, 9, false, false);
377  tmp[ 56] = ElmtConfig(eTriangle, 10, false, false);
378  tmp[ 57] = ElmtConfig(eQuadrilateral, 6, false, false);
379  tmp[ 58] = ElmtConfig(eQuadrilateral, 7, false, false);
380  tmp[ 59] = ElmtConfig(eQuadrilateral, 8, false, false);
381  tmp[ 60] = ElmtConfig(eQuadrilateral, 9, false, false);
382  tmp[ 61] = ElmtConfig(eQuadrilateral, 10, false, false);
383  tmp[ 62] = ElmtConfig(eSegment, 6, true, false);
384  tmp[ 63] = ElmtConfig(eSegment, 7, true, false);
385  tmp[ 64] = ElmtConfig(eSegment, 8, true, false);
386  tmp[ 65] = ElmtConfig(eSegment, 9, true, false);
387  tmp[ 66] = ElmtConfig(eSegment, 10, true, false);
388  tmp[ 71] = ElmtConfig(eTetrahedron, 6, true, true);
389  tmp[ 72] = ElmtConfig(eTetrahedron, 7, true, true);
390  tmp[ 73] = ElmtConfig(eTetrahedron, 8, true, true);
391  tmp[ 74] = ElmtConfig(eTetrahedron, 9, true, true);
392  tmp[ 75] = ElmtConfig(eTetrahedron, 10, true, true);
393  tmp[ 79] = ElmtConfig(eTetrahedron, 6, true, false);
394  tmp[ 80] = ElmtConfig(eTetrahedron, 7, true, false);
395  tmp[ 81] = ElmtConfig(eTetrahedron, 8, true, false);
396  tmp[ 82] = ElmtConfig(eTetrahedron, 9, true, false);
397  tmp[ 83] = ElmtConfig(eTetrahedron, 10, true, false);
398  tmp[ 90] = ElmtConfig(ePrism, 3, true, true);
399  tmp[ 91] = ElmtConfig(ePrism, 4, true, true);
400  tmp[ 92] = ElmtConfig(eHexahedron, 3, true, true);
401  tmp[ 93] = ElmtConfig(eHexahedron, 4, true, true);
402  tmp[ 94] = ElmtConfig(eHexahedron, 5, true, true);
403  tmp[ 95] = ElmtConfig(eHexahedron, 6, true, true);
404  tmp[ 96] = ElmtConfig(eHexahedron, 7, true, true);
405  tmp[ 97] = ElmtConfig(eHexahedron, 8, true, true);
406  tmp[ 98] = ElmtConfig(eHexahedron, 9, true, true);
407  tmp[ 99] = ElmtConfig(eHexahedron, 3, true, false);
408  tmp[100] = ElmtConfig(eHexahedron, 4, true, false);
409  tmp[101] = ElmtConfig(eHexahedron, 5, true, false);
410  tmp[102] = ElmtConfig(eHexahedron, 6, true, false);
411  tmp[103] = ElmtConfig(eHexahedron, 7, true, false);
412  tmp[104] = ElmtConfig(eHexahedron, 8, true, false);
413  tmp[105] = ElmtConfig(eHexahedron, 9, true, false);
414  tmp[106] = ElmtConfig(ePrism, 5, true, true);
415  tmp[107] = ElmtConfig(ePrism, 6, true, true);
416  tmp[108] = ElmtConfig(ePrism, 7, true, true);
417  tmp[109] = ElmtConfig(ePrism, 8, true, true);
418  tmp[110] = ElmtConfig(ePrism, 9, true, true);
419  tmp[111] = ElmtConfig(ePrism, 3, true, false);
420  tmp[112] = ElmtConfig(ePrism, 4, true, false);
421  tmp[113] = ElmtConfig(ePrism, 5, true, false);
422  tmp[114] = ElmtConfig(ePrism, 6, true, false);
423  tmp[115] = ElmtConfig(ePrism, 7, true, false);
424  tmp[116] = ElmtConfig(ePrism, 8, true, false);
425  tmp[117] = ElmtConfig(ePrism, 9, true, false);
426 
427  return tmp;
428  }
429  }
430 }
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2