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