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 
49 
50 #include "InputGmsh.h"
51 
52 using namespace Nektar::NekMeshUtils;
53 
54 namespace Nektar
55 {
56 namespace Utilities
57 {
58 
59 ModuleKey InputGmsh::className = GetModuleFactory().RegisterCreatorFunction(
60  ModuleKey(eInputModule, "msh"), InputGmsh::create, "Reads Gmsh msh file.");
61 
62 std::map<unsigned int, ElmtConfig> InputGmsh::elmMap = InputGmsh::GenElmMap();
63 
64 /**
65  * @brief Reorder a quadrilateral to appear in Nektar++ ordering from
66  * Gmsh.
67  */
68 std::vector<int> quadTensorNodeOrdering(const std::vector<int> &nodes, int n)
69 {
70  std::vector<int> nodeList;
71 
72  // Triangle
73  nodeList.resize(nodes.size());
74 
75  // Vertices and edges
76  nodeList[0] = nodes[0];
77  if (n > 1)
78  {
79  nodeList[n - 1] = nodes[1];
80  nodeList[n * n - 1] = nodes[2];
81  nodeList[n * (n - 1)] = nodes[3];
82  }
83  for (int i = 1; i < n - 1; ++i)
84  {
85  nodeList[i] = nodes[4 + i - 1];
86  }
87  for (int i = 1; i < n - 1; ++i)
88  {
89  nodeList[n * n - 1 - i] = nodes[4 + 2 * (n - 2) + i - 1];
90  }
91 
92  // Interior (recursion)
93  if (n > 2)
94  {
95  // Reorder interior nodes
96  std::vector<int> interior((n - 2) * (n - 2));
97  std::copy(
98  nodes.begin() + 4 + 4 * (n - 2), nodes.end(), interior.begin());
99  interior = quadTensorNodeOrdering(interior, n - 2);
100 
101  // Copy into full node list
102  for (int j = 1; j < n - 1; ++j)
103  {
104  nodeList[j * n] = nodes[4 + 3 * (n - 2) + n - 2 - j];
105  for (int i = 1; i < n - 1; ++i)
106  {
107  nodeList[j * n + i] = interior[(j - 1) * (n - 2) + (i - 1)];
108  }
109  nodeList[(j + 1) * n - 1] = nodes[4 + (n - 2) + j - 1];
110  }
111  }
112 
113  return nodeList;
114 }
115 
116 std::vector<int> triTensorNodeOrdering(const std::vector<int> &nodes, int n)
117 {
118  std::vector<int> nodeList;
119  int cnt2;
120 
121  nodeList.resize(nodes.size());
122 
123  // Vertices
124  nodeList[0] = nodes[0];
125  if (n > 1)
126  {
127  nodeList[n - 1] = nodes[1];
128  nodeList[n * (n + 1) / 2 - 1] = nodes[2];
129  }
130 
131  // Edges
132  int cnt = n;
133  for (int i = 1; i < n - 1; ++i)
134  {
135  nodeList[i] = nodes[3 + i - 1];
136  nodeList[cnt] = nodes[3 + 3 * (n - 2) - i];
137  nodeList[cnt + n - i - 1] = nodes[3 + (n - 2) + i - 1];
138  cnt += n - i;
139  }
140 
141  // Interior (recursion)
142  if (n > 3)
143  {
144  // Reorder interior nodes
145  std::vector<int> interior((n - 3) * (n - 2) / 2);
146  std::copy(
147  nodes.begin() + 3 + 3 * (n - 2), nodes.end(), interior.begin());
148  interior = triTensorNodeOrdering(interior, n - 3);
149 
150  // Copy into full node list
151  cnt = n;
152  cnt2 = 0;
153  for (int j = 1; j < n - 2; ++j)
154  {
155  for (int i = 0; i < n - j - 2; ++i)
156  {
157  nodeList[cnt + i + 1] = interior[cnt2 + i];
158  }
159  cnt += n - j;
160  cnt2 += n - 2 - j;
161  }
162  }
163 
164  return nodeList;
165 }
166 
167 /**
168  * @brief Set up InputGmsh object.
169  *
170  */
171 InputGmsh::InputGmsh(MeshSharedPtr m) : InputModule(m)
172 {
173 }
174 
176 {
177 }
178 
179 /**
180  * Gmsh file contains a list of nodes and their coordinates, along with
181  * a list of elements and those nodes which define them. We read in and
182  * store the list of nodes in #m_node and store the list of elements in
183  * #m_element. Each new element is supplied with a list of entries from
184  * #m_node which defines the element. Finally some mesh statistics are
185  * printed.
186  *
187  * @param pFilename Filename of Gmsh file to read.
188  */
190 {
191  // Open the file stream.
192  OpenStream();
193 
194  m_mesh->m_expDim = 0;
195  m_mesh->m_spaceDim = 0;
196  string line;
197  int nVertices = 0;
198  int nEntities = 0;
199  int elm_type = 0;
200  int prevId = -1;
201  int maxTagId = -1;
202 
204 
205  // This map takes each element ID and maps it to a permutation map
206  // that is required to take Gmsh element node orderings and map them
207  // to Nektar++ orderings.
208  boost::unordered_map<int, vector<int> > orderingMap;
209  boost::unordered_map<int, vector<int> >::iterator oIt;
210 
211  if (m_mesh->m_verbose)
212  {
213  cout << "InputGmsh: Start reading file..." << endl;
214  }
215 
216  while (!m_mshFile.eof())
217  {
218  getline(m_mshFile, line);
219  stringstream s(line);
220  string word;
221  s >> word;
222 
223  // Process nodes.
224  if (word == "$Nodes")
225  {
226  getline(m_mshFile, line);
227  stringstream s(line);
228  s >> nVertices;
229  int id = 0;
230  for (int i = 0; i < nVertices; ++i)
231  {
232  getline(m_mshFile, line);
233  stringstream st(line);
234  double x = 0, y = 0, z = 0;
235  st >> id >> x >> y >> z;
236 
237  if ((x * x) > 0.000001 && m_mesh->m_spaceDim < 1)
238  {
239  m_mesh->m_spaceDim = 1;
240  }
241  if ((y * y) > 0.000001 && m_mesh->m_spaceDim < 2)
242  {
243  m_mesh->m_spaceDim = 2;
244  }
245  if ((z * z) > 0.000001 && m_mesh->m_spaceDim < 3)
246  {
247  m_mesh->m_spaceDim = 3;
248  }
249 
250  id -= 1; // counter starts at 0
251 
252  if (id - prevId != 1)
253  {
254  cerr << "Gmsh vertex ids should be contiguous" << endl;
255  abort();
256  }
257  prevId = id;
258  m_mesh->m_node.push_back(
259  boost::shared_ptr<Node>(new Node(id, x, y, z)));
260  }
261  }
262  // Process elements
263  else if (word == "$Elements")
264  {
265  getline(m_mshFile, line);
266  stringstream s(line);
267  s >> nEntities;
268  for (int i = 0; i < nEntities; ++i)
269  {
270  getline(m_mshFile, line);
271  stringstream st(line);
272  int id = 0, num_tag = 0, num_nodes = 0;
273 
274  st >> id >> elm_type >> num_tag;
275  id -= 1; // counter starts at 0
276 
277  it = elmMap.find(elm_type);
278  if (it == elmMap.end())
279  {
280  cerr << "Error: element type " << elm_type
281  << " not supported" << endl;
282  abort();
283  }
284 
285  // Read element tags
286  vector<int> tags;
287  for (int j = 0; j < num_tag; ++j)
288  {
289  int tag = 0;
290  st >> tag;
291  tags.push_back(tag);
292  }
293  tags.resize(1);
294 
295  maxTagId = max(maxTagId, tags[0]);
296 
297  // Read element node list
298  vector<NodeSharedPtr> nodeList;
299  num_nodes = GetNnodes(elm_type);
300  for (int k = 0; k < num_nodes; ++k)
301  {
302  int node = 0;
303  st >> node;
304  node -= 1; // counter starts at 0
305  nodeList.push_back(m_mesh->m_node[node]);
306  }
307 
308  // Look up reordering.
309  oIt = orderingMap.find(elm_type);
310 
311  // If it's not created, then create it.
312  if (oIt == orderingMap.end())
313  {
314  oIt = orderingMap.insert(
315  make_pair(elm_type, CreateReordering(elm_type))).first;
316  }
317 
318  // Apply reordering map where necessary.
319  if (oIt->second.size() > 0)
320  {
321  vector<int> &mapping = oIt->second;
322  vector<NodeSharedPtr> tmp = nodeList;
323  for (int i = 0; i < mapping.size(); ++i)
324  {
325  nodeList[i] = tmp[mapping[i]];
326  }
327  }
328 
329  // Create element
331  it->second.m_e, it->second, nodeList, tags);
332 
333  // Determine mesh expansion dimension
334  if (E->GetDim() > m_mesh->m_expDim)
335  {
336  m_mesh->m_expDim = E->GetDim();
337  }
338  m_mesh->m_element[E->GetDim()].push_back(E);
339  }
340  }
341  }
342  m_mshFile.close();
343 
344  // Go through element and remap tags if necessary.
345  map<int, map<LibUtilities::ShapeType, int> > compMap;
346  map<int, map<LibUtilities::ShapeType, int> >::iterator cIt;
348 
349  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
350  {
351  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
352  LibUtilities::ShapeType type = el->GetConf().m_e;
353 
354  vector<int> tags = el->GetTagList();
355  int tag = tags[0];
356 
357  cIt = compMap.find(tag);
358 
359  if (cIt == compMap.end())
360  {
361  compMap[tag][type] = tag;
362  continue;
363  }
364 
365  // Reset tag for this element.
366  sIt = cIt->second.find(type);
367  if (sIt == cIt->second.end())
368  {
369  maxTagId++;
370  cIt->second[type] = maxTagId;
371  tags[0] = maxTagId;
372  el->SetTagList(tags);
373  }
374  else if (sIt->second != tag)
375  {
376  tags[0] = sIt->second;
377  el->SetTagList(tags);
378  }
379  }
380 
381  bool printInfo = false;
382  for (cIt = compMap.begin(); cIt != compMap.end(); ++cIt)
383  {
384  if (cIt->second.size() > 1)
385  {
386  printInfo = true;
387  break;
388  }
389  }
390 
391  if (printInfo)
392  {
393  cout << "Multiple elements in composite detected; remapped:" << endl;
394  for (cIt = compMap.begin(); cIt != compMap.end(); ++cIt)
395  {
396  if (cIt->second.size() > 1)
397  {
398  sIt = cIt->second.begin();
399  cout << "- Tag " << cIt->first << " => " << sIt->second << " ("
400  << LibUtilities::ShapeTypeMap[sIt->first] << ")";
401  sIt++;
402 
403  for (; sIt != cIt->second.end(); ++sIt)
404  {
405  cout << ", " << sIt->second << " ("
406  << LibUtilities::ShapeTypeMap[sIt->first] << ")";
407  }
408 
409  cout << endl;
410  }
411  }
412  }
413 
414  // Process rest of mesh.
415  ProcessVertices();
416  ProcessEdges();
417  ProcessFaces();
418  ProcessElements();
420 }
421 
422 /**
423  * For a given msh ID, return the corresponding number of nodes.
424  */
425 int InputGmsh::GetNnodes(unsigned int InputGmshEntity)
426 {
427  int nNodes;
429 
430  it = elmMap.find(InputGmshEntity);
431 
432  if (it == elmMap.end())
433  {
434  cerr << "Unknown element type " << InputGmshEntity << endl;
435  abort();
436  }
437 
438  switch (it->second.m_e)
439  {
441  nNodes = Point::GetNumNodes(it->second);
442  break;
444  nNodes = Line::GetNumNodes(it->second);
445  break;
447  nNodes = Triangle::GetNumNodes(it->second);
448  break;
450  nNodes = Quadrilateral::GetNumNodes(it->second);
451  break;
453  nNodes = Tetrahedron::GetNumNodes(it->second);
454  it->second.m_faceCurveType = LibUtilities::eNodalTriEvenlySpaced;
455  break;
457  nNodes = Pyramid::GetNumNodes(it->second);
458  break;
460  nNodes = Prism::GetNumNodes(it->second);
461  break;
463  nNodes = Hexahedron::GetNumNodes(it->second);
464  break;
465  default:
466  cerr << "Unknown element type!" << endl;
467  abort();
468  break;
469  }
470 
471  return nNodes;
472 }
473 
474 /**
475  * @brief Create a reordering map for a given element.
476  *
477  * Since Gmsh and Nektar++ have different vertex, edge and face
478  * orientations, we need to reorder the nodes in a Gmsh MSH file so that
479  * they work with the Nektar++ orderings, since this is what is used in
480  * the elements defined in the converter.
481  */
482 vector<int> InputGmsh::CreateReordering(unsigned int InputGmshEntity)
483 {
485 
486  it = elmMap.find(InputGmshEntity);
487 
488  if (it == elmMap.end())
489  {
490  cerr << "Unknown element type " << InputGmshEntity << endl;
491  abort();
492  }
493 
494  // For specific elements, call the appropriate function to perform
495  // the renumbering.
496  switch (it->second.m_e)
497  {
499  return TriReordering(it->second);
500  break;
502  return QuadReordering(it->second);
503  break;
505  return TetReordering(it->second);
506  break;
508  return PrismReordering(it->second);
509  break;
511  return HexReordering(it->second);
512  break;
513  default:
514  break;
515  }
516 
517  // Default: no reordering.
518  vector<int> returnVal;
519  return returnVal;
520 }
521 
522 /**
523  * @brief Create a reordering for triangles.
524  */
526 {
527  const int order = conf.m_order;
528  const int n = order - 1;
529 
530  // Copy vertices.
531  vector<int> mapping(3);
532  for (int i = 0; i < 3; ++i)
533  {
534  mapping[i] = i;
535  }
536 
537  if (order == 1)
538  {
539  return mapping;
540  }
541 
542  // Curvilinear edges.
543  mapping.resize(3 + 3 * n);
544 
545  for (int i = 3; i < 3 + 3 * n; ++i)
546  {
547  mapping[i] = i;
548  }
549 
550  if (!conf.m_faceNodes)
551  {
552  return mapping;
553  }
554 
555  // Interior nodes.
556  vector<int> interior(n * (n - 1) / 2);
557  for (int i = 0; i < interior.size(); ++i)
558  {
559  interior[i] = i + 3 + 3 * n;
560  }
561 
562  if (interior.size() > 0)
563  {
564  interior = triTensorNodeOrdering(interior, n - 1);
565  }
566 
567  mapping.insert(mapping.end(), interior.begin(), interior.end());
568  return mapping;
569 }
570 
571 /**
572  * @brief Create a reordering for quadrilaterals.
573  */
575 {
576  const int order = conf.m_order;
577  const int n = order - 1;
578 
579  // Copy vertices.
580  vector<int> mapping(4);
581  for (int i = 0; i < 4; ++i)
582  {
583  mapping[i] = i;
584  }
585 
586  if (order == 1)
587  {
588  return mapping;
589  }
590 
591  // Curvilinear edges.
592  mapping.resize(4 + 4 * n);
593 
594  for (int i = 4; i < 4 + 4 * n; ++i)
595  {
596  mapping[i] = i;
597  }
598 
599  if (!conf.m_faceNodes)
600  {
601  return mapping;
602  }
603 
604  // Interior nodes.
605  vector<int> interior(n * n);
606  for (int i = 0; i < interior.size(); ++i)
607  {
608  interior[i] = i + 4 + 4 * n;
609  }
610 
611  if (interior.size() > 0)
612  {
613  interior = quadTensorNodeOrdering(interior, n);
614  }
615  mapping.insert(mapping.end(), interior.begin(), interior.end());
616  return mapping;
617 }
618 
619 /**
620  * @brief Create a reordering for tetrahedra.
621  */
623 {
624  const int order = conf.m_order;
625  const int n = order - 1;
626  const int n2 = n * (n - 1) / 2;
627 
628  int i, j;
629  vector<int> mapping(4);
630 
631  // Copy vertices.
632  for (i = 0; i < 4; ++i)
633  {
634  mapping[i] = i;
635  }
636 
637  if (order == 1)
638  {
639  return mapping;
640  }
641 
642  // Curvilinear edges.
643  mapping.resize(4 + 6 * n);
644 
645  // Curvilinear edges.
646  static int gmshToNekEdge[6] = {0, 1, 2, 3, 5, 4};
647  static int gmshToNekRev[6] = {0, 0, 1, 1, 1, 1};
648 
649  // Reorder edges.
650  int offset, cnt = 4;
651  for (i = 0; i < 6; ++i)
652  {
653  offset = 4 + n * gmshToNekEdge[i];
654 
655  if (gmshToNekRev[i])
656  {
657  for (int j = 0; j < n; ++j)
658  {
659  mapping[offset + n - j - 1] = cnt++;
660  }
661  }
662  else
663  {
664  for (int j = 0; j < n; ++j)
665  {
666  mapping[offset + j] = cnt++;
667  }
668  }
669  }
670 
671  if (conf.m_faceNodes == false || n2 == 0)
672  {
673  return mapping;
674  }
675 
676  // Curvilinear faces.
677  mapping.resize(4 + 6 * n + 4 * n2);
678 
679  static int gmshToNekFace[4] = {0, 1, 3, 2};
680 
681  vector<int> triVertId(3);
682  triVertId[0] = 0;
683  triVertId[1] = 1;
684  triVertId[2] = 2;
685 
686  // Loop over Gmsh faces
687  for (i = 0; i < 4; ++i)
688  {
689  int face = gmshToNekFace[i];
690  int offset2 = 4 + 6 * n + i * n2;
691  offset = 4 + 6 * n + face * n2;
692 
693  // Create a list of interior face nodes for this face only.
694  vector<int> faceNodes(n2);
695  vector<int> toAlign(3);
696  for (j = 0; j < n2; ++j)
697  {
698  faceNodes[j] = offset2 + j;
699  }
700 
701  // Now get the reordering of this face, which puts Gmsh
702  // recursive ordering into Nektar++ row-by-row order.
703  vector<int> tmp = triTensorNodeOrdering(faceNodes, n - 1);
704  HOTriangle<int> hoTri(triVertId, tmp);
705 
706  // Apply reorientation
707  if (i == 0 || i == 2)
708  {
709  // Triangle verts {0,2,1} --> {0,1,2}
710  toAlign[0] = 0;
711  toAlign[1] = 2;
712  toAlign[2] = 1;
713  hoTri.Align(toAlign);
714  }
715  else if (i == 3)
716  {
717  // Triangle verts {1,2,0} --> {0,1,2}
718  toAlign[0] = 1;
719  toAlign[1] = 2;
720  toAlign[2] = 0;
721  hoTri.Align(toAlign);
722  }
723 
724  // Fill in mapping.
725  for (j = 0; j < n2; ++j)
726  {
727  mapping[offset + j] = hoTri.surfVerts[j];
728  }
729  }
730 
731  return mapping;
732 }
733 
734 /**
735  * @brief Create a reordering for prisms.
736  *
737  * Note that whilst Gmsh MSH files have the capability to support
738  * high-order prisms, presently Gmsh does not seem to be capable of
739  * generating higher than second-order prismatic meshes, so most of the
740  * following is untested.
741  */
743 {
744  const int order = conf.m_order;
745  const int n = order - 1;
746 
747  int i;
748  vector<int> mapping(6);
749 
750  // To get from Gmsh -> Nektar++ prism, coordinates axes are
751  // different; need to mirror in the triangular faces, and then
752  // reorder vertices to make ordering anticlockwise on base quad.
753  static int gmshToNekVerts[6] = {3, 4, 1, 0, 5, 2};
754 
755  for (i = 0; i < 6; ++i)
756  {
757  mapping[i] = gmshToNekVerts[i];
758  }
759 
760  if (order == 1)
761  {
762  return mapping;
763  }
764 
765  // Curvilinear edges.
766  mapping.resize(6 + 9 * n);
767 
768  static int gmshToNekEdge[9] = {2, 7, 3, 6, 1, 8, 0, 4, 5};
769  static int gmshToNekRev[9] = {1, 0, 1, 0, 1, 1, 0, 0, 0};
770 
771  // Reorder edges.
772  int offset, cnt = 6;
773  for (i = 0; i < 9; ++i)
774  {
775  offset = 6 + n * gmshToNekEdge[i];
776 
777  if (gmshToNekRev[i])
778  {
779  for (int j = 0; j < n; ++j)
780  {
781  mapping[offset + n - j - 1] = cnt++;
782  }
783  }
784  else
785  {
786  for (int j = 0; j < n; ++j)
787  {
788  mapping[offset + j] = cnt++;
789  }
790  }
791  }
792 
793  if (conf.m_faceNodes == false)
794  {
795  return mapping;
796  }
797 
798  if (order > 2)
799  {
800  cerr << "Gmsh prisms of order > 2 with face curvature "
801  << "not supported in NekMesh (or indeed Gmsh at"
802  << "time of writing)." << endl;
803  abort();
804  }
805 
806  mapping.resize(18);
807  mapping[15] = 15;
808  mapping[16] = 17;
809  mapping[17] = 16;
810 
811  return mapping;
812 }
813 
814 /**
815  * @brief Create a reordering for hexahedra.
816  */
818 {
819  const int order = conf.m_order;
820  const int n = order - 1;
821  const int n2 = n * n;
822  int i, j, k;
823 
824  vector<int> mapping;
825 
826  // Map taking Gmsh edges to Nektar++ edges.
827  static int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10};
828 
829  // Push back vertices.
830  mapping.resize(8);
831  for (i = 0; i < 8; ++i)
832  {
833  mapping[i] = i;
834  }
835 
836  if (order == 1)
837  {
838  return mapping;
839  }
840 
841  // Curvilinear edges
842  mapping.resize(8 + 12 * n);
843 
844  // Reorder edges.
845  int cnt = 8, offset;
846  for (i = 0; i < 12; ++i)
847  {
848  int edge = gmshToNekEdge[i];
849  offset = 8 + n * abs(edge);
850 
851  if (edge < 0)
852  {
853  for (int j = 0; j < n; ++j)
854  {
855  mapping[offset + n - j - 1] = cnt++;
856  }
857  }
858  else
859  {
860  for (int j = 0; j < n; ++j)
861  {
862  mapping[offset + j] = cnt++;
863  }
864  }
865  }
866 
867  if (conf.m_faceNodes == false || n2 == 0)
868  {
869  return mapping;
870  }
871 
872  // Curvilinear face nodes.
873  mapping.resize(8 + 12 * n + 6 * n2);
874 
875  // Map which takes Gmsh -> Nektar++ faces in the local element.
876  static int gmsh2NekFace[6] = {0, 1, 4, 2, 3, 5};
877 
878  // Map which defines orientation between Gmsh and Nektar++ faces.
879  StdRegions::Orientation faceOrient[6] = {
885  StdRegions::eDir1FwdDir1_Dir2FwdDir2};
886 
887  for (i = 0; i < 6; ++i)
888  {
889  int face = gmsh2NekFace[i];
890  int offset2 = 8 + 12 * n + i * n2;
891  offset = 8 + 12 * n + face * n2;
892 
893  // Create a list of interior face nodes for this face only.
894  vector<int> faceNodes(n2);
895  for (j = 0; j < n2; ++j)
896  {
897  faceNodes[j] = offset2 + j;
898  }
899 
900  // Now get the reordering of this face, which puts Gmsh
901  // recursive ordering into Nektar++ row-by-row order.
902  vector<int> tmp = quadTensorNodeOrdering(faceNodes, n);
903 
904  // Finally reorient the face according to the geometry
905  // differences.
906  if (faceOrient[i] == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
907  {
908  // Orientation is the same, just copy.
909  for (j = 0; j < n2; ++j)
910  {
911  mapping[offset + j] = tmp[j];
912  }
913  }
914  else if (faceOrient[i] == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
915  {
916  // Tranposed faces
917  for (j = 0; j < n; ++j)
918  {
919  for (k = 0; k < n; ++k)
920  {
921  mapping[offset + j * n + k] = tmp[k * n + j];
922  }
923  }
924  }
925  else if (faceOrient[i] == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
926  {
927  for (j = 0; j < n; ++j)
928  {
929  for (k = 0; k < n; ++k)
930  {
931  mapping[offset + j * n + k] = tmp[j * n + (n - k - 1)];
932  }
933  }
934  }
935  }
936 
937  if (conf.m_volumeNodes == false)
938  {
939  return mapping;
940  }
941 
942  const int totPoints = (order + 1) * (order + 1) * (order + 1);
943  mapping.resize(totPoints);
944 
945  // TODO: Fix ordering of volume nodes.
946  for (i = 8 + 12 * n + 6 * n2; i < totPoints; ++i)
947  {
948  mapping[i] = i;
949  }
950 
951  return mapping;
952 }
953 
954 /*
955  * @brief Populate the element map #elmMap.
956  *
957  * This function primarily populates the element mapping #elmMap,
958  * which takes a msh ID used by Gmsh and translates to element type,
959  * element order and whether the element is incomplete (i.e. whether
960  * it contains solely boundary nodes, or just face nodes). Note that
961  * some of these elements, such as prisms of order >= 3, cannot yet be
962  * generated by Gmsh.
963  */
964 std::map<unsigned int, ElmtConfig> InputGmsh::GenElmMap()
965 {
966  using namespace LibUtilities;
967  std::map<unsigned int, ElmtConfig> tmp;
968 
969  // Elmt type, order, face, volume
970  tmp[ 1] = ElmtConfig(eSegment, 1, true, true);
971  tmp[ 2] = ElmtConfig(eTriangle, 1, true, true);
972  tmp[ 3] = ElmtConfig(eQuadrilateral, 1, true, true);
973  tmp[ 4] = ElmtConfig(eTetrahedron, 1, true, true);
974  tmp[ 5] = ElmtConfig(eHexahedron, 1, true, true);
975  tmp[ 6] = ElmtConfig(ePrism, 1, true, true);
976  tmp[ 7] = ElmtConfig(ePyramid, 1, true, true);
977  tmp[ 8] = ElmtConfig(eSegment, 2, true, true);
978  tmp[ 9] = ElmtConfig(eTriangle, 2, true, true);
979  tmp[ 10] = ElmtConfig(eQuadrilateral, 2, true, true);
980  tmp[ 11] = ElmtConfig(eTetrahedron, 2, true, true);
981  tmp[ 12] = ElmtConfig(eHexahedron, 2, true, true);
982  tmp[ 13] = ElmtConfig(ePrism, 2, true, true);
983  tmp[ 14] = ElmtConfig(ePyramid, 2, true, true);
984  tmp[ 15] = ElmtConfig(ePoint, 1, true, false);
985  tmp[ 16] = ElmtConfig(eQuadrilateral, 2, false, false);
986  tmp[ 17] = ElmtConfig(eHexahedron, 2, false, false);
987  tmp[ 18] = ElmtConfig(ePrism, 2, false, false);
988  tmp[ 19] = ElmtConfig(ePyramid, 2, false, false);
989  tmp[ 20] = ElmtConfig(eTriangle, 3, false, false);
990  tmp[ 21] = ElmtConfig(eTriangle, 3, true, false);
991  tmp[ 22] = ElmtConfig(eTriangle, 4, false, false);
992  tmp[ 23] = ElmtConfig(eTriangle, 4, true, false);
993  tmp[ 24] = ElmtConfig(eTriangle, 5, false, false);
994  tmp[ 25] = ElmtConfig(eTriangle, 5, true, false);
995  tmp[ 26] = ElmtConfig(eSegment, 3, true, false);
996  tmp[ 27] = ElmtConfig(eSegment, 4, true, false);
997  tmp[ 28] = ElmtConfig(eSegment, 5, true, false);
998  tmp[ 29] = ElmtConfig(eTetrahedron, 3, true, true);
999  tmp[ 30] = ElmtConfig(eTetrahedron, 4, true, true);
1000  tmp[ 31] = ElmtConfig(eTetrahedron, 5, true, true);
1001  tmp[ 32] = ElmtConfig(eTetrahedron, 4, true, false);
1002  tmp[ 33] = ElmtConfig(eTetrahedron, 5, true, false);
1003  tmp[ 36] = ElmtConfig(eQuadrilateral, 3, true, false);
1004  tmp[ 37] = ElmtConfig(eQuadrilateral, 4, true, false);
1005  tmp[ 38] = ElmtConfig(eQuadrilateral, 5, true, false);
1006  tmp[ 39] = ElmtConfig(eQuadrilateral, 3, false, false);
1007  tmp[ 40] = ElmtConfig(eQuadrilateral, 4, false, false);
1008  tmp[ 41] = ElmtConfig(eQuadrilateral, 5, false, false);
1009  tmp[ 42] = ElmtConfig(eTriangle, 6, true, false);
1010  tmp[ 43] = ElmtConfig(eTriangle, 7, true, false);
1011  tmp[ 44] = ElmtConfig(eTriangle, 8, true, false);
1012  tmp[ 45] = ElmtConfig(eTriangle, 9, true, false);
1013  tmp[ 46] = ElmtConfig(eTriangle, 10, true, false);
1014  tmp[ 47] = ElmtConfig(eQuadrilateral, 6, true, false);
1015  tmp[ 48] = ElmtConfig(eQuadrilateral, 7, true, false);
1016  tmp[ 49] = ElmtConfig(eQuadrilateral, 8, true, false);
1017  tmp[ 50] = ElmtConfig(eQuadrilateral, 9, true, false);
1018  tmp[ 51] = ElmtConfig(eQuadrilateral, 10, true, false);
1019  tmp[ 52] = ElmtConfig(eTriangle, 6, false, false);
1020  tmp[ 53] = ElmtConfig(eTriangle, 7, false, false);
1021  tmp[ 54] = ElmtConfig(eTriangle, 8, false, false);
1022  tmp[ 55] = ElmtConfig(eTriangle, 9, false, false);
1023  tmp[ 56] = ElmtConfig(eTriangle, 10, false, false);
1024  tmp[ 57] = ElmtConfig(eQuadrilateral, 6, false, false);
1025  tmp[ 58] = ElmtConfig(eQuadrilateral, 7, false, false);
1026  tmp[ 59] = ElmtConfig(eQuadrilateral, 8, false, false);
1027  tmp[ 60] = ElmtConfig(eQuadrilateral, 9, false, false);
1028  tmp[ 61] = ElmtConfig(eQuadrilateral, 10, false, false);
1029  tmp[ 62] = ElmtConfig(eSegment, 6, true, false);
1030  tmp[ 63] = ElmtConfig(eSegment, 7, true, false);
1031  tmp[ 64] = ElmtConfig(eSegment, 8, true, false);
1032  tmp[ 65] = ElmtConfig(eSegment, 9, true, false);
1033  tmp[ 66] = ElmtConfig(eSegment, 10, true, false);
1034  tmp[ 71] = ElmtConfig(eTetrahedron, 6, true, true);
1035  tmp[ 72] = ElmtConfig(eTetrahedron, 7, true, true);
1036  tmp[ 73] = ElmtConfig(eTetrahedron, 8, true, true);
1037  tmp[ 74] = ElmtConfig(eTetrahedron, 9, true, true);
1038  tmp[ 75] = ElmtConfig(eTetrahedron, 10, true, true);
1039  tmp[ 79] = ElmtConfig(eTetrahedron, 6, true, false);
1040  tmp[ 80] = ElmtConfig(eTetrahedron, 7, true, false);
1041  tmp[ 81] = ElmtConfig(eTetrahedron, 8, true, false);
1042  tmp[ 82] = ElmtConfig(eTetrahedron, 9, true, false);
1043  tmp[ 83] = ElmtConfig(eTetrahedron, 10, true, false);
1044  tmp[ 90] = ElmtConfig(ePrism, 3, true, true);
1045  tmp[ 91] = ElmtConfig(ePrism, 4, true, true);
1046  tmp[ 92] = ElmtConfig(eHexahedron, 3, true, true);
1047  tmp[ 93] = ElmtConfig(eHexahedron, 4, true, true);
1048  tmp[ 94] = ElmtConfig(eHexahedron, 5, true, true);
1049  tmp[ 95] = ElmtConfig(eHexahedron, 6, true, true);
1050  tmp[ 96] = ElmtConfig(eHexahedron, 7, true, true);
1051  tmp[ 97] = ElmtConfig(eHexahedron, 8, true, true);
1052  tmp[ 98] = ElmtConfig(eHexahedron, 9, true, true);
1053  tmp[ 99] = ElmtConfig(eHexahedron, 3, true, false);
1054  tmp[100] = ElmtConfig(eHexahedron, 4, true, false);
1055  tmp[101] = ElmtConfig(eHexahedron, 5, true, false);
1056  tmp[102] = ElmtConfig(eHexahedron, 6, true, false);
1057  tmp[103] = ElmtConfig(eHexahedron, 7, true, false);
1058  tmp[104] = ElmtConfig(eHexahedron, 8, true, false);
1059  tmp[105] = ElmtConfig(eHexahedron, 9, true, false);
1060  tmp[106] = ElmtConfig(ePrism, 5, true, true);
1061  tmp[107] = ElmtConfig(ePrism, 6, true, true);
1062  tmp[108] = ElmtConfig(ePrism, 7, true, true);
1063  tmp[109] = ElmtConfig(ePrism, 8, true, true);
1064  tmp[110] = ElmtConfig(ePrism, 9, true, true);
1065  tmp[111] = ElmtConfig(ePrism, 3, true, false);
1066  tmp[112] = ElmtConfig(ePrism, 4, true, false);
1067  tmp[113] = ElmtConfig(ePrism, 5, true, false);
1068  tmp[114] = ElmtConfig(ePrism, 6, true, false);
1069  tmp[115] = ElmtConfig(ePrism, 7, true, false);
1070  tmp[116] = ElmtConfig(ePrism, 8, true, false);
1071  tmp[117] = ElmtConfig(ePrism, 9, true, false);
1072  tmp[118] = ElmtConfig(ePyramid, 3, true, true);
1073  tmp[119] = ElmtConfig(ePyramid, 4, true, true);
1074  tmp[120] = ElmtConfig(ePyramid, 5, true, true);
1075  tmp[121] = ElmtConfig(ePyramid, 6, true, true);
1076  tmp[122] = ElmtConfig(ePyramid, 7, true, true);
1077  tmp[123] = ElmtConfig(ePyramid, 8, true, true);
1078  tmp[124] = ElmtConfig(ePyramid, 9, true, true);
1079  tmp[125] = ElmtConfig(ePyramid, 3, true, false);
1080  tmp[126] = ElmtConfig(ePyramid, 4, true, false);
1081  tmp[127] = ElmtConfig(ePyramid, 5, true, false);
1082  tmp[128] = ElmtConfig(ePyramid, 6, true, false);
1083  tmp[129] = ElmtConfig(ePyramid, 7, true, false);
1084  tmp[130] = ElmtConfig(ePyramid, 7, true, false);
1085  tmp[131] = ElmtConfig(ePyramid, 8, true, false);
1086 
1087  return tmp;
1088 }
1089 
1090 }
1091 }
bool m_faceNodes
Denotes whether the element contains face nodes. For 2D elements, if this is true then the element co...
Definition: Element.h:89
std::vector< int > triTensorNodeOrdering(const std::vector< int > &nodes, int n)
Definition: InputGmsh.cpp:116
Basic information about an element.
Definition: Element.h:58
std::ifstream m_mshFile
Input stream.
pair< ModuleType, string > ModuleKey
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
void Align(std::vector< int > vertId)
Align this surface to a given vertex ID.
Definition: Triangle.h:138
std::vector< int > HexReordering(ElmtConfig conf)
Create a reordering for hexahedra.
Definition: InputGmsh.cpp:817
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a quadrilateral.
std::vector< int > quadTensorNodeOrdering(const std::vector< int > &nodes, int n)
Reorder a quadrilateral to appear in Nektar++ ordering from Gmsh.
Definition: InputGmsh.cpp:68
std::vector< int > TriReordering(ElmtConfig conf)
Create a reordering for triangles.
Definition: InputGmsh.cpp:525
STL namespace.
std::vector< T > surfVerts
The triangle surface vertices – templated so that this can either be nodes or IDs.
Definition: Triangle.h:68
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
Represents a point in the domain.
Definition: Node.h:60
const char *const ShapeTypeMap[]
Definition: ShapeType.hpp:66
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
unsigned int m_order
Order of the element.
Definition: Element.h:96
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a line.
Definition: Line.cpp:115
virtual void ProcessVertices()
Extract element vertices.
std::vector< int > QuadReordering(ElmtConfig conf)
Create a reordering for quadrilaterals.
Definition: InputGmsh.cpp:574
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a prism.
Definition: Prism.cpp:218
std::vector< int > CreateReordering(unsigned int InputGmshEntity)
Create a reordering map for a given element.
Definition: InputGmsh.cpp:482
bool m_volumeNodes
Denotes whether the element contains volume (i.e. interior) nodes. These are not supported by either ...
Definition: Element.h:94
virtual void ProcessElements()
Generate element IDs.
std::vector< int > PrismReordering(ElmtConfig conf)
Create a reordering for prisms.
Definition: InputGmsh.cpp:742
virtual void ProcessComposites()
Generate composites.
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a triangle.
Definition: Triangle.cpp:144
A lightweight struct for dealing with high-order triangle alignment.
Definition: Triangle.h:53
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a hexahedron.
Definition: Hexahedron.cpp:182
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::vector< int > TetReordering(ElmtConfig conf)
Create a reordering for tetrahedra.
Definition: InputGmsh.cpp:622
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:137
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a tetrahedron.
int GetNnodes(unsigned int InputGmshEntity)
Definition: InputGmsh.cpp:425
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a pyramid.
Definition: Pyramid.cpp:175
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a point (i.e. return 1).
Definition: Point.cpp:66
2D Evenly-spaced points on a Triangle
Definition: PointsType.h:70
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
static std::map< unsigned int, ElmtConfig > GenElmMap()
Definition: InputGmsh.cpp:964
void OpenStream()
Open a file for input.
ModuleFactory & GetModuleFactory()
Abstract base class for input modules.
static std::map< unsigned int, ElmtConfig > elmMap
Definition: InputGmsh.h:68
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215