Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 <boost/tuple/tuple.hpp>
41 
51 
52 #include "InputGmsh.h"
53 
54 using namespace Nektar::NekMeshUtils;
55 
56 namespace Nektar
57 {
58 namespace Utilities
59 {
60 
61 ModuleKey InputGmsh::className = GetModuleFactory().RegisterCreatorFunction(
62  ModuleKey(eInputModule, "msh"), InputGmsh::create, "Reads Gmsh msh file.");
63 
64 std::map<unsigned int, ElmtConfig> InputGmsh::elmMap = InputGmsh::GenElmMap();
65 
66 /**
67  * @brief Reorder Gmsh nodes so that they appear in a tensor product format
68  * suitable for the interior of a Nektar++ quadrilateral.
69  *
70  * For an example, consider a second order quadrilateral. This routine will
71  * produce a permutation map that applies the following permutation:
72  *
73  * ```
74  * 3---6---2 6---7---8
75  * | | | |
76  * 7 8 5 ---> 3 4 5
77  * | | | |
78  * 0---4---1 0---1---2
79  *
80  * Gmsh tensor-product
81  * ```
82  *
83  * We assume that Gmsh uses a recursive ordering system, so that interior nodes
84  * are reordered by calling this function recursively.
85  *
86  * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format.
87  * @param n The number of nodes in one coordinate direction. If this is
88  * zero we assume no reordering needs to be done and return the
89  * identity permutation.
90  *
91  * @return The nodes vector in tensor-product ordering.
92  */
93 std::vector<int> quadTensorNodeOrdering(const std::vector<int> &nodes, int n)
94 {
95  std::vector<int> nodeList;
96 
97  // Triangle
98  nodeList.resize(nodes.size());
99 
100  // Vertices and edges
101  nodeList[0] = nodes[0];
102  if (n > 1)
103  {
104  nodeList[n - 1] = nodes[1];
105  nodeList[n * n - 1] = nodes[2];
106  nodeList[n * (n - 1)] = nodes[3];
107  }
108  for (int i = 1; i < n - 1; ++i)
109  {
110  nodeList[i] = nodes[4 + i - 1];
111  }
112  for (int i = 1; i < n - 1; ++i)
113  {
114  nodeList[n * n - 1 - i] = nodes[4 + 2 * (n - 2) + i - 1];
115  }
116 
117  // Interior (recursion)
118  if (n > 2)
119  {
120  // Reorder interior nodes
121  std::vector<int> interior((n - 2) * (n - 2));
122  std::copy(
123  nodes.begin() + 4 + 4 * (n - 2), nodes.end(), interior.begin());
124  interior = quadTensorNodeOrdering(interior, n - 2);
125 
126  // Copy into full node list
127  for (int j = 1; j < n - 1; ++j)
128  {
129  nodeList[j * n] = nodes[4 + 3 * (n - 2) + n - 2 - j];
130  for (int i = 1; i < n - 1; ++i)
131  {
132  nodeList[j * n + i] = interior[(j - 1) * (n - 2) + (i - 1)];
133  }
134  nodeList[(j + 1) * n - 1] = nodes[4 + (n - 2) + j - 1];
135  }
136  }
137 
138  return nodeList;
139 }
140 
141 /**
142  * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
143  * suitable for the interior of a Nektar++ triangle.
144  *
145  * For an example, consider a third-order triangle. This routine will produce a
146  * permutation map that applies the following permutation:
147  *
148  * ```
149  * 2 9
150  * | \ | \
151  * 7 6 7 7
152  * | \ | \
153  * 8 9 5 4 5 6
154  * | \ | \
155  * 0---3---4---1 0---1---2---3
156  *
157  * Gmsh tensor-product
158  * ```
159  *
160  * We assume that Gmsh uses a recursive ordering system, so that interior nodes
161  * are reordered by calling this function recursively.
162  *
163  * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format.
164  * @param n The number of nodes in one coordinate direction. If this is
165  * zero we assume no reordering needs to be done and return the
166  * identity permutation.
167  *
168  * @return The nodes vector in tensor-product ordering.
169  */
170 std::vector<int> triTensorNodeOrdering(const std::vector<int> &nodes, int n)
171 {
172  std::vector<int> nodeList;
173  int cnt2;
174 
175  nodeList.resize(nodes.size());
176 
177  // Vertices
178  nodeList[0] = nodes[0];
179  if (n > 1)
180  {
181  nodeList[n - 1] = nodes[1];
182  nodeList[n * (n + 1) / 2 - 1] = nodes[2];
183  }
184 
185  // Edges
186  int cnt = n;
187  for (int i = 1; i < n - 1; ++i)
188  {
189  nodeList[i] = nodes[3 + i - 1];
190  nodeList[cnt] = nodes[3 + 3 * (n - 2) - i];
191  nodeList[cnt + n - i - 1] = nodes[3 + (n - 2) + i - 1];
192  cnt += n - i;
193  }
194 
195  // Interior (recursion)
196  if (n > 3)
197  {
198  // Reorder interior nodes
199  std::vector<int> interior((n - 3) * (n - 2) / 2);
200  std::copy(
201  nodes.begin() + 3 + 3 * (n - 2), nodes.end(), interior.begin());
202  interior = triTensorNodeOrdering(interior, n - 3);
203 
204  // Copy into full node list
205  cnt = n;
206  cnt2 = 0;
207  for (int j = 1; j < n - 2; ++j)
208  {
209  for (int i = 0; i < n - j - 2; ++i)
210  {
211  nodeList[cnt + i + 1] = interior[cnt2 + i];
212  }
213  cnt += n - j;
214  cnt2 += n - 2 - j;
215  }
216  }
217 
218  return nodeList;
219 }
220 
221 typedef boost::tuple<int, int, int> Mode;
222 struct cmpop
223 {
224  bool operator()(Mode const &a, Mode const &b) const
225  {
226  if (a.get<0>() < b.get<0>())
227  {
228  return true;
229  }
230  if (a.get<0>() > b.get<0>())
231  {
232  return false;
233  }
234  if (a.get<1>() < b.get<1>())
235  {
236  return true;
237  }
238  if (a.get<1>() > b.get<1>())
239  {
240  return false;
241  }
242  if (a.get<2>() < b.get<2>())
243  {
244  return true;
245  }
246 
247  return false;
248  }
249 };
250 
251 /**
252  * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
253  * suitable for the interior of a Nektar++ tetrahedron.
254  *
255  * For an example, consider a second order tetrahedron. This routine will
256  * produce a permutation map that applies the following permutation:
257  *
258  * ```
259  * 2 5
260  * ,/|`\ ,/|`\
261  * ,/ | `\ ,/ | `\
262  * ,6 '. `5 ,3 '. `4
263  * ,/ 8 `\ ,/ 8 `\
264  * ,/ | `\ ,/ | `\
265  * 0--------4--'.--------1 0--------1--'.--------2
266  * `\. | ,/ `\. | ,/
267  * `\. | ,9 `\. | ,7
268  * `7. '. ,/ `6. '. ,/
269  * `\. |/ `\. |/
270  * `3 `9
271  *
272  * Gmsh tensor-product
273  * ```
274  *
275  * We assume that Gmsh uses a recursive ordering system, so that interior nodes
276  * are reordered by calling this function recursively.
277  *
278  * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format.
279  * @param n The number of nodes in one coordinate direction. If this is
280  * zero we assume no reordering needs to be done and return the
281  * identity permutation.
282  *
283  * @return The nodes vector in tensor-product ordering.
284  */
285 std::vector<int> tetTensorNodeOrdering(const std::vector<int> &nodes, int n)
286 {
287  std::vector<int> nodeList;
288  int nTri = n*(n+1)/2;
289  int nTet = n*(n+1)*(n+2)/6;
290 
291  nodeList.resize(nodes.size());
292  nodeList[0] = nodes[0];
293 
294  if (n == 1)
295  {
296  return nodeList;
297  }
298 
299  // Vertices
300  nodeList[n - 1] = nodes[1];
301  nodeList[nTri - 1] = nodes[2];
302  nodeList[nTet - 1] = nodes[3];
303 
304  if (n == 2)
305  {
306  return nodeList;
307  }
308 
309  // Set up a map that takes (a,b,c) -> m to help us figure out where things
310  // are inside the tetrahedron.
311  std::map<Mode, int, cmpop> tmp;
312 
313  for (int k = 0, cnt = 0; k < n; ++k)
314  {
315  for (int j = 0; j < n - k; ++j)
316  {
317  for (int i = 0; i < n - k - j; ++i)
318  {
319  tmp[Mode(i,j,k)] = cnt++;
320  }
321  }
322  }
323 
324  // Edges
325  for (int i = 1; i < n-1; ++i)
326  {
327  int eI = i-1;
328  nodeList[tmp[Mode(i,0,0)]] = nodes[4 + eI];
329  nodeList[tmp[Mode(n-1-i,i,0)]] = nodes[4 + (n-2) + eI];
330  nodeList[tmp[Mode(0,n-1-i,0)]] = nodes[4 + 2*(n-2) + eI];
331  nodeList[tmp[Mode(0,0,n-1-i)]] = nodes[4 + 3*(n-2) + eI];
332  nodeList[tmp[Mode(0,i,n-1-i)]] = nodes[4 + 4*(n-2) + eI];
333  nodeList[tmp[Mode(i,0,n-1-i)]] = nodes[4 + 5*(n-2) + eI];
334  }
335 
336  if (n == 3)
337  {
338  return nodeList;
339  }
340 
341  // For faces, we use the triTensorNodeOrdering routine to make our lives
342  // slightly easier.
343  int nFacePts = (n-3)*(n-2)/2;
344 
345  // Grab face points and reorder into a tensor-product type format
346  vector<vector<int> > tmpNodes(4);
347  int offset = 4 + 6*(n-2);
348 
349  for (int i = 0; i < 4; ++i)
350  {
351  tmpNodes[i].resize(nFacePts);
352  for (int j = 0; j < nFacePts; ++j)
353  {
354  tmpNodes[i][j] = nodes[offset++];
355  }
356  tmpNodes[i] = triTensorNodeOrdering(tmpNodes[i], n-3);
357  }
358 
359  if (n > 4)
360  {
361  // Now align faces
362  vector<int> triVertId(3), toAlign(3);
363  triVertId[0] = 0;
364  triVertId[1] = 1;
365  triVertId[2] = 2;
366 
367  // Faces 0,2: triangle verts {0,2,1} --> {0,1,2}
368  HOTriangle<int> hoTri(triVertId, tmpNodes[0]);
369  toAlign[0] = 0;
370  toAlign[1] = 2;
371  toAlign[2] = 1;
372 
373  hoTri.Align(toAlign);
374  tmpNodes[0] = hoTri.surfVerts;
375 
376  hoTri.surfVerts = tmpNodes[2];
377  hoTri.Align(toAlign);
378  tmpNodes[2] = hoTri.surfVerts;
379 
380  // Face 3: triangle verts {1,2,0} --> {0,1,2}
381  toAlign[0] = 1;
382  toAlign[1] = 2;
383  toAlign[2] = 0;
384 
385  hoTri.surfVerts = tmpNodes[3];
386  hoTri.Align(toAlign);
387  tmpNodes[3] = hoTri.surfVerts;
388  }
389 
390  // Now apply faces. Note that faces 3 and 2 are swapped between Gmsh and
391  // Nektar++ order.
392  for (int j = 1, cnt = 0; j < n-2; ++j)
393  {
394  for (int i = 1; i < n-j-1; ++i, ++cnt)
395  {
396  nodeList[tmp[Mode(i,j,0)]] = tmpNodes[0][cnt];
397  nodeList[tmp[Mode(i,0,j)]] = tmpNodes[1][cnt];
398  nodeList[tmp[Mode(n-1-i-j,i,j)]] = tmpNodes[3][cnt];
399  nodeList[tmp[Mode(0,i,j)]] = tmpNodes[2][cnt];
400  }
401  }
402 
403  if (n == 4)
404  {
405  return nodeList;
406  }
407 
408  // Finally, recurse on interior volume
409  vector<int> intNodes, tmpInt;
410  for (int i = offset; i < nTet; ++i)
411  {
412  intNodes.push_back(nodes[i]);
413  }
414  tmpInt = tetTensorNodeOrdering(intNodes, n-4);
415 
416  for (int k = 1, cnt = 0; k < n - 2; ++k)
417  {
418  for (int j = 1; j < n - k - 1; ++j)
419  {
420  for (int i = 1; i < n - k - j - 1; ++i)
421  {
422  nodeList[tmp[Mode(i,j,k)]] = tmpInt[cnt++];
423  }
424  }
425  }
426 
427  return nodeList;
428 }
429 
430 /**
431  * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
432  * suitable for the interior of a Nektar++ prism. This routine is specifically
433  * designed for *interior* Gmsh nodes only.
434  *
435  * Prisms are a bit of a special case, in that interior nodes are heterogeneous
436  * in the number of points they use. As an example, for a second-order interior
437  * of a fourth-order prism, this routine calculates the mapping taking us
438  *
439  * ```
440  * ,6 ,8
441  * .-' | \ .-' | \
442  * ,-3 | \ ,-5 | \
443  * ,-' | \ | \ ,-' | \ | \
444  * 2 | \7-------8 2 | \6-------7
445  * | \ | ,-' \ ,-' | \ | ,-' \ ,-'
446  * | \,5------,0' | \,3------,4'
447  * |,-' \ ,-' |,-' \ ,-'
448  * 1-------4' 0-------1'
449  *
450  * Gmsh tensor-product
451  * ```
452  *
453  * @todo Make this routine work for higher orders. The Gmsh ordering seems
454  * a little different than other elements, therefore the functionality is
455  * (for now) hard-coded to orders less than 5.
456  *
457  * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format.
458  * @param n The number of nodes in one coordinate direction. If this is
459  * zero we assume no reordering needs to be done and return the
460  * identity permutation.
461  *
462  * @return The nodes vector in tensor-product ordering.
463  */
464 std::vector<int> prismTensorNodeOrdering(const std::vector<int> &nodes, int n)
465 {
466  std::vector<int> nodeList;
467 
468  if (n == 0)
469  {
470  return nodeList;
471  }
472 
473  nodeList.resize(nodes.size());
474 
475  if (n == 2)
476  {
477  nodeList[0] = nodes[1];
478  nodeList[1] = nodes[0];
479  return nodeList;
480  }
481 
482  // For some reason, this ordering is different. Whereas Gmsh usually orders
483  // vertices first, followed by edges, this ordering goes VVE-VVE-VVE for the
484  // three edges that contain edge-interior information, but also vertex
485  // orientation is not the same as the original Gmsh prism. The below is done
486  // by looking at the ordering by hand and is hence why we don't support
487  // order > 4 right now.
488  if (n == 3)
489  {
490  nodeList[0] = nodes[1];
491  nodeList[1] = nodes[4];
492  nodeList[2] = nodes[2];
493  nodeList[3] = nodes[5];
494  nodeList[4] = nodes[0];
495  nodeList[5] = nodes[3];
496  nodeList[6] = nodes[7];
497  nodeList[7] = nodes[8];
498  nodeList[8] = nodes[6];
499  }
500 
501  ASSERTL0(n < 4, "Prism Gmsh input and output is incomplete for orders "
502  "larger than 4");
503 
504  return nodeList;
505 }
506 
507 /**
508  * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
509  * suitable for the interior of a Nektar++ hexahedron.
510  *
511  * For an example, consider a second order hexahedron. This routine will produce
512  * a permutation map that applies the following permutation:
513  *
514  * ```
515  * 3----13----2 6----7-----8
516  * |\ |\ |\ |\
517  * |15 24 | 14 |15 16 | 18
518  * 9 \ 20 11 \ 3 \ 4 5 \
519  * | 7----19+---6 | 24---25+---26
520  * |22 | 26 | 23| |12 | 13 | 14|
521  * 0---+-8----1 | 0---+-1----2 |
522  * \ 17 25 \ 18 \ 21 22 \ 23
523  * 10 | 21 12| 9 | 10 11|
524  * \| \| \| \|
525  * 4----16----5 18---19----20
526  *
527  * Gmsh tensor-product
528  * ```
529  *
530  * We assume that Gmsh uses a recursive ordering system, so that interior nodes
531  * are reordered by calling this function recursively.
532  *
533  * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format.
534  * @param n The number of nodes in one coordinate direction. If this is
535  * zero we assume no reordering needs to be done and return the
536  * identity permutation.
537  *
538  * @return The nodes vector in tensor-product ordering.
539  */
540 std::vector<int> hexTensorNodeOrdering(const std::vector<int> &nodes, int n)
541 {
542  int i, j, k;
543  std::vector<int> nodeList;
544 
545  nodeList.resize(nodes.size());
546  nodeList[0] = nodes[0];
547 
548  if (n == 1)
549  {
550  return nodeList;
551  }
552 
553  // Vertices: same order as Nektar++
554  nodeList[n - 1] = nodes[1];
555  nodeList[n*n -1] = nodes[2];
556  nodeList[n*(n-1)] = nodes[3];
557  nodeList[n*n*(n-1)] = nodes[4];
558  nodeList[n - 1 + n*n*(n-1)] = nodes[5];
559  nodeList[n*n -1 + n*n*(n-1)] = nodes[6];
560  nodeList[n*(n-1) + n*n*(n-1)] = nodes[7];
561 
562  if (n == 2)
563  {
564  return nodeList;
565  }
566 
567  int hexEdges[12][2] = {
568  { 0, 1 }, { n-1, n }, { n*n-1, -1 }, { n*(n-1), -n },
569  { 0, n*n }, { n-1, n*n }, { n*n - 1, n*n }, { n*(n-1), n*n },
570  { n*n*(n-1), 1 }, { n*n*(n-1) + n-1, n }, { n*n*n-1, -1 },
571  { n*n*(n-1) + n*(n-1), -n }
572  };
573  int hexFaces[6][3] = {
574  { 0, 1, n }, { 0, 1, n*n }, { n-1, n, n*n },
575  { n*(n-1), 1, n*n }, { 0, n, n*n }, { n*n*(n-1), 1, n }
576  };
577  int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10};
578 
579  // Edges
580  int offset = 8;
581  for (int i = 0; i < 12; ++i)
582  {
583  int e = abs(gmshToNekEdge[i]);
584 
585  if (gmshToNekEdge[i] >= 0)
586  {
587  for (int j = 1; j < n-1; ++j)
588  {
589  nodeList[hexEdges[e][0] + j*hexEdges[e][1]] = nodes[offset++];
590  }
591  }
592  else
593  {
594  for (int j = 1; j < n-1; ++j)
595  {
596  nodeList[hexEdges[e][0] + (n-j-1)*hexEdges[e][1]] = nodes[offset++];
597  }
598  }
599  }
600 
601  // Faces
602  int gmsh2NekFace[6] = {0, 1, 4, 2, 3, 5};
603 
604  // Map which defines orientation between Gmsh and Nektar++ faces.
605  StdRegions::Orientation faceOrient[6] = {
611  StdRegions::eDir1FwdDir1_Dir2FwdDir2};
612 
613  for (i = 0; i < 6; ++i)
614  {
615  int n2 = (n-2)*(n-2);
616  int face = gmsh2NekFace[i];
617  offset = 8 + 12 * (n-2) + i * n2;
618 
619  // Create a list of interior face nodes for this face only.
620  vector<int> faceNodes(n2);
621  for (j = 0; j < n2; ++j)
622  {
623  faceNodes[j] = nodes[offset + j];
624  }
625 
626  // Now get the reordering of this face, which puts Gmsh
627  // recursive ordering into Nektar++ row-by-row order.
628  faceNodes = quadTensorNodeOrdering(faceNodes, n-2);
629  vector<int> tmp(n2);
630 
631  // Finally reorient the face according to the geometry
632  // differences.
633  if (faceOrient[i] == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
634  {
635  // Orientation is the same, just copy.
636  tmp = faceNodes;
637  }
638  else if (faceOrient[i] == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
639  {
640  // Tranposed faces
641  for (j = 0; j < n-2; ++j)
642  {
643  for (k = 0; k < n-2; ++k)
644  {
645  tmp[j * (n-2) + k] = faceNodes[k * (n-2) + j];
646  }
647  }
648  }
649  else if (faceOrient[i] == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
650  {
651  for (j = 0; j < n-2; ++j)
652  {
653  for (k = 0; k < n-2; ++k)
654  {
655  tmp[j * (n-2) + k] = faceNodes[j * (n-2) + (n - k - 3)];
656  }
657  }
658  }
659 
660  // Now put this into the right place in the output array
661  for (k = 1; k < n-1; ++k)
662  {
663  for (j = 1; j < n-1; ++j)
664  {
665  nodeList[hexFaces[face][0] + j*hexFaces[face][1] + k*hexFaces[face][2]]
666  = faceNodes[(k-1)*(n-2) + j-1];
667  }
668  }
669  }
670 
671  // Finally, recurse on interior volume
672  vector<int> intNodes, tmpInt;
673  for (int i = 8 + 12 * (n-2) + 6 * (n-2) * (n-2); i < n*n*n; ++i)
674  {
675  intNodes.push_back(nodes[i]);
676  }
677 
678  if (intNodes.size())
679  {
680  tmpInt = hexTensorNodeOrdering(intNodes, n-2);
681  for (int k = 1, cnt = 0; k < n - 1; ++k)
682  {
683  for (int j = 1; j < n - 1; ++j)
684  {
685  for (int i = 1; i < n - 1; ++i)
686  {
687  nodeList[i + j * n + k * n * n] = tmpInt[cnt++];
688  }
689  }
690  }
691  }
692 
693  return nodeList;
694 }
695 
696 
697 /**
698  * @brief Set up InputGmsh object.
699  *
700  */
701 InputGmsh::InputGmsh(MeshSharedPtr m) : InputModule(m)
702 {
703 }
704 
706 {
707 }
708 
709 /**
710  * Gmsh file contains a list of nodes and their coordinates, along with a list
711  * of elements and those nodes which define them. We read in and store the list
712  * of nodes in #m_node and store the list of elements in #m_element. Each new
713  * element is supplied with a list of entries from m_node which defines the
714  * element. Finally some mesh statistics are printed.
715  *
716  * @param pFilename Filename of Gmsh file to read.
717  */
719 {
720  // Open the file stream.
721  OpenStream();
722 
723  m_mesh->m_expDim = 0;
724  m_mesh->m_spaceDim = 0;
725  string line;
726  int nVertices = 0;
727  int nEntities = 0;
728  int elm_type = 0;
729  int prevId = -1;
730  int maxTagId = -1;
731 
733 
734  // This map takes each element ID and maps it to a permutation map
735  // that is required to take Gmsh element node orderings and map them
736  // to Nektar++ orderings.
737  boost::unordered_map<int, vector<int> > orderingMap;
738  boost::unordered_map<int, vector<int> >::iterator oIt;
739 
740  if (m_mesh->m_verbose)
741  {
742  cout << "InputGmsh: Start reading file..." << endl;
743  }
744 
745  while (!m_mshFile.eof())
746  {
747  getline(m_mshFile, line);
748  stringstream s(line);
749  string word;
750  s >> word;
751 
752  // Process nodes.
753  if (word == "$Nodes")
754  {
755  getline(m_mshFile, line);
756  stringstream s(line);
757  s >> nVertices;
758  int id = 0;
759  for (int i = 0; i < nVertices; ++i)
760  {
761  getline(m_mshFile, line);
762  stringstream st(line);
763  double x = 0, y = 0, z = 0;
764  st >> id >> x >> y >> z;
765 
766  if ((x * x) > 0.000001 && m_mesh->m_spaceDim < 1)
767  {
768  m_mesh->m_spaceDim = 1;
769  }
770  if ((y * y) > 0.000001 && m_mesh->m_spaceDim < 2)
771  {
772  m_mesh->m_spaceDim = 2;
773  }
774  if ((z * z) > 0.000001 && m_mesh->m_spaceDim < 3)
775  {
776  m_mesh->m_spaceDim = 3;
777  }
778 
779  id -= 1; // counter starts at 0
780 
781  if (id - prevId != 1)
782  {
783  cerr << "Gmsh vertex ids should be contiguous" << endl;
784  abort();
785  }
786  prevId = id;
787  m_mesh->m_node.push_back(
788  boost::shared_ptr<Node>(new Node(id, x, y, z)));
789  }
790  }
791  // Process elements
792  else if (word == "$Elements")
793  {
794  getline(m_mshFile, line);
795  stringstream s(line);
796  s >> nEntities;
797  for (int i = 0; i < nEntities; ++i)
798  {
799  getline(m_mshFile, line);
800  stringstream st(line);
801  int id = 0, num_tag = 0, num_nodes = 0;
802 
803  st >> id >> elm_type >> num_tag;
804  id -= 1; // counter starts at 0
805 
806  it = elmMap.find(elm_type);
807  if (it == elmMap.end())
808  {
809  cerr << "Error: element type " << elm_type
810  << " not supported" << endl;
811  abort();
812  }
813 
814  // Read element tags
815  vector<int> tags;
816  for (int j = 0; j < num_tag; ++j)
817  {
818  int tag = 0;
819  st >> tag;
820  tags.push_back(tag);
821  }
822  tags.resize(1);
823 
824  maxTagId = max(maxTagId, tags[0]);
825 
826  // Read element node list
827  vector<NodeSharedPtr> nodeList;
828  num_nodes = GetNnodes(elm_type);
829  for (int k = 0; k < num_nodes; ++k)
830  {
831  int node = 0;
832  st >> node;
833  node -= 1; // counter starts at 0
834  nodeList.push_back(m_mesh->m_node[node]);
835  }
836 
837  // Look up reordering.
838  oIt = orderingMap.find(elm_type);
839 
840  // If it's not created, then create it.
841  if (oIt == orderingMap.end())
842  {
843  oIt = orderingMap.insert(
844  make_pair(elm_type, CreateReordering(elm_type))).first;
845  }
846 
847  // Apply reordering map where necessary.
848  if (oIt->second.size() > 0)
849  {
850  vector<int> &mapping = oIt->second;
851  vector<NodeSharedPtr> tmp = nodeList;
852  for (int i = 0; i < mapping.size(); ++i)
853  {
854  nodeList[i] = tmp[mapping[i]];
855  }
856  }
857 
858  // Create element
860  it->second.m_e, it->second, nodeList, tags);
861 
862  // Determine mesh expansion dimension
863  if (E->GetDim() > m_mesh->m_expDim)
864  {
865  m_mesh->m_expDim = E->GetDim();
866  }
867  m_mesh->m_element[E->GetDim()].push_back(E);
868  }
869  }
870  }
871  m_mshFile.reset();
872 
873  // Go through element and remap tags if necessary.
874  map<int, map<LibUtilities::ShapeType, int> > compMap;
875  map<int, map<LibUtilities::ShapeType, int> >::iterator cIt;
877 
878  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
879  {
880  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
881  LibUtilities::ShapeType type = el->GetConf().m_e;
882 
883  vector<int> tags = el->GetTagList();
884  int tag = tags[0];
885 
886  cIt = compMap.find(tag);
887 
888  if (cIt == compMap.end())
889  {
890  compMap[tag][type] = tag;
891  continue;
892  }
893 
894  // Reset tag for this element.
895  sIt = cIt->second.find(type);
896  if (sIt == cIt->second.end())
897  {
898  maxTagId++;
899  cIt->second[type] = maxTagId;
900  tags[0] = maxTagId;
901  el->SetTagList(tags);
902  }
903  else if (sIt->second != tag)
904  {
905  tags[0] = sIt->second;
906  el->SetTagList(tags);
907  }
908  }
909 
910  bool printInfo = false;
911  for (cIt = compMap.begin(); cIt != compMap.end(); ++cIt)
912  {
913  if (cIt->second.size() > 1)
914  {
915  printInfo = true;
916  break;
917  }
918  }
919 
920  if (printInfo)
921  {
922  cout << "Multiple elements in composite detected; remapped:" << endl;
923  for (cIt = compMap.begin(); cIt != compMap.end(); ++cIt)
924  {
925  if (cIt->second.size() > 1)
926  {
927  sIt = cIt->second.begin();
928  cout << "- Tag " << cIt->first << " => " << sIt->second << " ("
929  << LibUtilities::ShapeTypeMap[sIt->first] << ")";
930  sIt++;
931 
932  for (; sIt != cIt->second.end(); ++sIt)
933  {
934  cout << ", " << sIt->second << " ("
935  << LibUtilities::ShapeTypeMap[sIt->first] << ")";
936  }
937 
938  cout << endl;
939  }
940  }
941  }
942 
943  // Process rest of mesh.
944  ProcessVertices();
945  ProcessEdges();
946  ProcessFaces();
947  ProcessElements();
949 }
950 
951 /**
952  * For a given msh ID, return the corresponding number of nodes.
953  */
954 int InputGmsh::GetNnodes(unsigned int InputGmshEntity)
955 {
956  int nNodes;
958 
959  it = elmMap.find(InputGmshEntity);
960 
961  if (it == elmMap.end())
962  {
963  cerr << "Unknown element type " << InputGmshEntity << endl;
964  abort();
965  }
966 
967  switch (it->second.m_e)
968  {
970  nNodes = Point::GetNumNodes(it->second);
971  break;
973  nNodes = Line::GetNumNodes(it->second);
974  break;
976  nNodes = Triangle::GetNumNodes(it->second);
977  break;
979  nNodes = Quadrilateral::GetNumNodes(it->second);
980  break;
982  nNodes = Tetrahedron::GetNumNodes(it->second);
983  it->second.m_faceCurveType = LibUtilities::eNodalTriEvenlySpaced;
984  break;
986  nNodes = Pyramid::GetNumNodes(it->second);
987  break;
989  nNodes = Prism::GetNumNodes(it->second);
990  break;
992  nNodes = Hexahedron::GetNumNodes(it->second);
993  break;
994  default:
995  cerr << "Unknown element type!" << endl;
996  abort();
997  break;
998  }
999 
1000  return nNodes;
1001 }
1002 
1003 /**
1004  * @brief Create a reordering map for a given element.
1005  *
1006  * Since Gmsh and Nektar++ have different vertex, edge and face
1007  * orientations, we need to reorder the nodes in a Gmsh MSH file so that
1008  * they work with the Nektar++ orderings, since this is what is used in
1009  * the elements defined in the converter.
1010  */
1011 vector<int> InputGmsh::CreateReordering(unsigned int InputGmshEntity)
1012 {
1014 
1015  it = elmMap.find(InputGmshEntity);
1016 
1017  if (it == elmMap.end())
1018  {
1019  cerr << "Unknown element type " << InputGmshEntity << endl;
1020  abort();
1021  }
1022 
1023  // For specific elements, call the appropriate function to perform
1024  // the renumbering.
1025  switch (it->second.m_e)
1026  {
1028  return TriReordering(it->second);
1029  break;
1031  return QuadReordering(it->second);
1032  break;
1034  return TetReordering(it->second);
1035  break;
1036  case LibUtilities::ePrism:
1037  return PrismReordering(it->second);
1038  break;
1040  return HexReordering(it->second);
1041  break;
1043  return LineReordering(it->second);
1044  default:
1045  break;
1046  }
1047 
1048  // Default: no reordering.
1049  vector<int> returnVal;
1050  return returnVal;
1051 }
1052 
1054 {
1055  const int order = conf.m_order;
1056 
1057  vector<int> mapping(order+1);
1058  for(int i = 0; i < order+1; i++)
1059  {
1060  mapping[i] = i;
1061  }
1062 
1063  return mapping;
1064 }
1065 
1066 /**
1067  * @brief Create a reordering for triangles.
1068  */
1070 {
1071  const int order = conf.m_order;
1072  const int n = order - 1;
1073 
1074  // Copy vertices.
1075  vector<int> mapping(3);
1076  for (int i = 0; i < 3; ++i)
1077  {
1078  mapping[i] = i;
1079  }
1080 
1081  if (order == 1)
1082  {
1083  return mapping;
1084  }
1085 
1086  // Curvilinear edges.
1087  mapping.resize(3 + 3 * n);
1088 
1089  for (int i = 3; i < 3 + 3 * n; ++i)
1090  {
1091  mapping[i] = i;
1092  }
1093 
1094  if (!conf.m_faceNodes)
1095  {
1096  return mapping;
1097  }
1098 
1099  // Interior nodes.
1100  vector<int> interior(n * (n - 1) / 2);
1101  for (int i = 0; i < interior.size(); ++i)
1102  {
1103  interior[i] = i + 3 + 3 * n;
1104  }
1105 
1106  if (interior.size() > 0)
1107  {
1108  interior = triTensorNodeOrdering(interior, n - 1);
1109  }
1110 
1111  mapping.insert(mapping.end(), interior.begin(), interior.end());
1112  return mapping;
1113 }
1114 
1115 /**
1116  * @brief Create a reordering for quadrilaterals.
1117  */
1119 {
1120  const int order = conf.m_order;
1121  const int n = order - 1;
1122 
1123  // Copy vertices.
1124  vector<int> mapping(4);
1125  for (int i = 0; i < 4; ++i)
1126  {
1127  mapping[i] = i;
1128  }
1129 
1130  if (order == 1)
1131  {
1132  return mapping;
1133  }
1134 
1135  // Curvilinear edges.
1136  mapping.resize(4 + 4 * n);
1137 
1138  for (int i = 4; i < 4 + 4 * n; ++i)
1139  {
1140  mapping[i] = i;
1141  }
1142 
1143  if (!conf.m_faceNodes)
1144  {
1145  return mapping;
1146  }
1147 
1148  // Interior nodes.
1149  vector<int> interior(n * n);
1150  for (int i = 0; i < interior.size(); ++i)
1151  {
1152  interior[i] = i + 4 + 4 * n;
1153  }
1154 
1155  if (interior.size() > 0)
1156  {
1157  interior = quadTensorNodeOrdering(interior, n);
1158  }
1159  mapping.insert(mapping.end(), interior.begin(), interior.end());
1160  return mapping;
1161 }
1162 
1163 /**
1164  * @brief Create a reordering for tetrahedra.
1165  */
1167 {
1168  const int order = conf.m_order;
1169  const int n = order - 1;
1170  const int n2 = n * (n - 1) / 2;
1171 
1172  int i, j;
1173  vector<int> mapping(4);
1174 
1175  // Copy vertices.
1176  for (i = 0; i < 4; ++i)
1177  {
1178  mapping[i] = i;
1179  }
1180 
1181  if (order == 1)
1182  {
1183  return mapping;
1184  }
1185 
1186  // Curvilinear edges.
1187  mapping.resize(4 + 6 * n);
1188 
1189  // Curvilinear edges.
1190  static int gmshToNekEdge[6] = {0, 1, 2, 3, 5, 4};
1191  static int gmshToNekRev[6] = {0, 0, 1, 1, 1, 1};
1192 
1193  // Reorder edges.
1194  int offset, cnt = 4;
1195  for (i = 0; i < 6; ++i)
1196  {
1197  offset = 4 + n * gmshToNekEdge[i];
1198 
1199  if (gmshToNekRev[i])
1200  {
1201  for (int j = 0; j < n; ++j)
1202  {
1203  mapping[offset + n - j - 1] = cnt++;
1204  }
1205  }
1206  else
1207  {
1208  for (int j = 0; j < n; ++j)
1209  {
1210  mapping[offset + j] = cnt++;
1211  }
1212  }
1213  }
1214 
1215  if (conf.m_faceNodes == false || n2 == 0)
1216  {
1217  return mapping;
1218  }
1219 
1220  // Curvilinear faces.
1221  mapping.resize(4 + 6 * n + 4 * n2);
1222 
1223  static int gmshToNekFace[4] = {0, 1, 3, 2};
1224 
1225  vector<int> triVertId(3);
1226  triVertId[0] = 0;
1227  triVertId[1] = 1;
1228  triVertId[2] = 2;
1229 
1230  // Loop over Gmsh faces
1231  for (i = 0; i < 4; ++i)
1232  {
1233  int face = gmshToNekFace[i];
1234  int offset2 = 4 + 6 * n + i * n2;
1235  offset = 4 + 6 * n + face * n2;
1236 
1237  // Create a list of interior face nodes for this face only.
1238  vector<int> faceNodes(n2);
1239  vector<int> toAlign(3);
1240  for (j = 0; j < n2; ++j)
1241  {
1242  faceNodes[j] = offset2 + j;
1243  }
1244 
1245  // Now get the reordering of this face, which puts Gmsh
1246  // recursive ordering into Nektar++ row-by-row order.
1247  vector<int> tmp = triTensorNodeOrdering(faceNodes, n - 1);
1248  HOTriangle<int> hoTri(triVertId, tmp);
1249 
1250  // Apply reorientation
1251  if (i == 0 || i == 2)
1252  {
1253  // Triangle verts {0,2,1} --> {0,1,2}
1254  toAlign[0] = 0;
1255  toAlign[1] = 2;
1256  toAlign[2] = 1;
1257  hoTri.Align(toAlign);
1258  }
1259  else if (i == 3)
1260  {
1261  // Triangle verts {1,2,0} --> {0,1,2}
1262  toAlign[0] = 1;
1263  toAlign[1] = 2;
1264  toAlign[2] = 0;
1265  hoTri.Align(toAlign);
1266  }
1267 
1268  // Fill in mapping.
1269  for (j = 0; j < n2; ++j)
1270  {
1271  mapping[offset + j] = hoTri.surfVerts[j];
1272  }
1273  }
1274 
1275  if (conf.m_volumeNodes == false)
1276  {
1277  return mapping;
1278  }
1279 
1280  const int nInt = (order - 3) * (order - 2) * (order - 1) / 6;
1281  if (nInt <= 0)
1282  {
1283  return mapping;
1284  }
1285 
1286  if (nInt == 1)
1287  {
1288  mapping.push_back(mapping.size());
1289  return mapping;
1290  }
1291 
1292  int ntot = (order+1)*(order+2)*(order+3)/6;
1293  vector<int> interior;
1294 
1295  for (int i = 4 + 6 * n + 4 * n2; i < ntot; ++i)
1296  {
1297  interior.push_back(i);
1298  }
1299 
1300  if (interior.size() > 0)
1301  {
1302  interior = tetTensorNodeOrdering(interior, order-3);
1303  }
1304 
1305  mapping.insert(mapping.end(), interior.begin(), interior.end());
1306 
1307  return mapping;
1308 }
1309 
1310 /**
1311  * @brief Create a reordering for prisms.
1312  *
1313  * Note that whilst Gmsh MSH files have the capability to support
1314  * high-order prisms, presently Gmsh does not seem to be capable of
1315  * generating higher than second-order prismatic meshes, so most of the
1316  * following is untested.
1317  */
1319 {
1320  const int order = conf.m_order;
1321  const int n = order - 1;
1322 
1323  int i;
1324  vector<int> mapping(6);
1325 
1326  // To get from Gmsh -> Nektar++ prism, coordinates axes are
1327  // different; need to mirror in the triangular faces, and then
1328  // reorder vertices to make ordering anticlockwise on base quad.
1329  static int nekToGmshVerts[6] = {3, 4, 1, 0, 5, 2};
1330  // Inverse (gmsh vert -> Nektar++ vert): 3->0, 4->1, 1->2, 0->3, 5->4, 2->5
1331 
1332  for (i = 0; i < 6; ++i)
1333  {
1334  mapping[i] = nekToGmshVerts[i];
1335  }
1336 
1337  if (order == 1)
1338  {
1339  return mapping;
1340  }
1341 
1342  // Curvilinear edges.
1343  mapping.resize(6 + 9 * n);
1344 
1345  // Nektar++ edges:
1346  // 0 = 1 = 2 = 3 = 4 = 5 = 6 = 7 = 8 =
1347  // {0,1}, {1,2}, {3,2}, {0,3}, {0,4}, {1,4}, {2,5}, {3,5}, {4,5}
1348  // Gmsh edges (from Geo/MPrism.h of Gmsh source):
1349  // {0,1}, {0,2}, {0,3}, {1,2}, {1,4}, {2,5}, {3,4}, {3,5}, {4,5}
1350  // Apply inverse of gmshToNekVerts map:
1351  // {3,2}, {3,5}, {3,0}, {2,5}, {2,1}, {5,4}, {0,1}, {0,4}, {1,4}
1352  // Nektar++ mapping (negative indicates reverse orientation):
1353  // = 2 = 7 = -3 = 6 = -1 = -8 = 0 = 4 = 5
1354  static int gmshToNekEdge[9] = {2, 7, 3, 6, 1, 8, 0, 4, 5};
1355  static int gmshToNekRev[9] = {0, 0, 1, 0, 1, 1, 0, 0, 0};
1356 
1357  // Reorder edges.
1358  int offset, cnt = 6;
1359  for (i = 0; i < 9; ++i)
1360  {
1361  offset = 6 + n * gmshToNekEdge[i];
1362 
1363  if (gmshToNekRev[i])
1364  {
1365  for (int j = 0; j < n; ++j)
1366  {
1367  mapping[offset + n - j - 1] = cnt++;
1368  }
1369  }
1370  else
1371  {
1372  for (int j = 0; j < n; ++j)
1373  {
1374  mapping[offset + j] = cnt++;
1375  }
1376  }
1377  }
1378 
1379  if (conf.m_faceNodes == false)
1380  {
1381  return mapping;
1382  }
1383 
1384  int nTriInt = n * (n - 1) / 2;
1385  int nQuadInt = n * n;
1386 
1387  // Nektar++ faces:
1388  // 0 = {0,1,2,3}, 1 = {0,1,4}, 2 = {1,2,5,4}, 3 = {3,2,5}, 4 = {0,3,5,4}
1389  // Gmsh faces (from Geo/MPrism.h of Gmsh source):
1390  // {0,2,1}, {3,4,5}, {0,1,4,3}, {0,3,5,2}, {1,2,5,4}
1391  // Apply inverse of gmshToNekVerts map:
1392  // {3,5,2}, {0,1,4}, {3,2,1,0}, {3,0,4,5}, {2,5,4,1}
1393  // = 3 = 1 = 0 = 4 = 2
1394  // This gives gmsh -> Nektar++ faces:
1395  static int gmshToNekFace[5] = {3, 1, 0, 4, 2};
1396 
1397  // Face offsets
1398  vector<int> offsets(5), offsets2(5);
1399  offset = 6 + 9 * n;
1400 
1401  // Offsets in the gmsh order: face ordering is TTQQQ
1402  offsets[0] = offset;
1403  offsets[1] = offset + nTriInt;
1404  offsets[2] = offset + 2 * nTriInt;
1405  offsets[3] = offset + 2 * nTriInt + nQuadInt;
1406  offsets[4] = offset + 2 * nTriInt + 2 * nQuadInt;
1407 
1408  // Offsets in the Nektar++ order: face ordering is QTQTQ
1409  offsets2[0] = offset;
1410  offsets2[1] = offset + nQuadInt;
1411  offsets2[2] = offset + nQuadInt + nTriInt;
1412  offsets2[3] = offset + 2 * nQuadInt + nTriInt;
1413  offsets2[4] = offset + 2 * nQuadInt + 2 * nTriInt;
1414 
1415  mapping.resize(6 + 9 * n + 3 * nQuadInt + 2 * nTriInt);
1416 
1417  offset = 6 + 9 * n;
1418 
1419  vector<int> triVertId(3);
1420  triVertId[0] = 0;
1421  triVertId[1] = 1;
1422  triVertId[2] = 2;
1423 
1424  vector<int> quadVertId(4);
1425  quadVertId[0] = 0;
1426  quadVertId[1] = 1;
1427  quadVertId[2] = 2;
1428  quadVertId[3] = 3;
1429 
1430  for (i = 0; i < 5; ++i)
1431  {
1432  int face = gmshToNekFace[i];
1433  int offset2 = offsets[i];
1434  offset = offsets2[face];
1435 
1436  bool tri = i < 2;
1437  int nFacePts = tri ? nTriInt : nQuadInt;
1438 
1439  if (nFacePts == 0)
1440  {
1441  continue;
1442  }
1443 
1444  // Create a list of interior face nodes for this face only.
1445  vector<int> faceNodes(nFacePts);
1446  vector<int> toAlign(tri ? 3 : 4);
1447  for (int j = 0; j < nFacePts; ++j)
1448  {
1449  faceNodes[j] = offset2 + j;
1450  }
1451 
1452  if (tri)
1453  {
1454  // Now get the reordering of this face, which puts Gmsh
1455  // recursive ordering into Nektar++ row-by-row order.
1456  vector<int> tmp = triTensorNodeOrdering(faceNodes, n - 1);
1457  HOTriangle<int> hoTri(triVertId, tmp);
1458 
1459  // Apply reorientation
1460  if (i == 0)
1461  {
1462  // Triangle verts {0,2,1} --> {0,1,2}
1463  toAlign[0] = 0;
1464  toAlign[1] = 2;
1465  toAlign[2] = 1;
1466  hoTri.Align(toAlign);
1467  }
1468 
1469  // Fill in mapping.
1470  for (int j = 0; j < nTriInt; ++j)
1471  {
1472  mapping[offset + j] = hoTri.surfVerts[j];
1473  }
1474  }
1475  else
1476  {
1477  vector<int> tmp = quadTensorNodeOrdering(faceNodes, n);
1478  HOQuadrilateral<int> hoQuad(quadVertId, tmp);
1479 
1480  // Apply reorientation
1481  if (i == 2)
1482  {
1483  toAlign[0] = 3;
1484  toAlign[1] = 2;
1485  toAlign[2] = 1;
1486  toAlign[3] = 0;
1487  }
1488  else if (i == 3)
1489  {
1490  toAlign[0] = 1;
1491  toAlign[1] = 0;
1492  toAlign[2] = 3;
1493  toAlign[3] = 2;
1494  }
1495  else if (i == 4)
1496  {
1497  toAlign[0] = 3;
1498  toAlign[1] = 0;
1499  toAlign[2] = 1;
1500  toAlign[3] = 2;
1501  }
1502 
1503  hoQuad.Align(toAlign);
1504 
1505  // Fill in mapping.
1506  for (int j = 0; j < nQuadInt; ++j)
1507  {
1508  mapping[offset + j] = hoQuad.surfVerts[j];
1509  }
1510  }
1511  }
1512 
1513  if (conf.m_volumeNodes == false)
1514  {
1515  return mapping;
1516  }
1517 
1518  // Interior points
1519  offset = offsets[4] + nQuadInt;
1520  vector<int> intPoints, tmp;
1521 
1522  for (int i = offset; i < (order+1) * (order+1) * (order+2) / 2; ++i)
1523  {
1524  intPoints.push_back(i);
1525  }
1526 
1527  // Reorder interior points
1528  tmp = prismTensorNodeOrdering(intPoints, order - 1);
1529  mapping.insert(mapping.end(), tmp.begin(), tmp.end());
1530 
1531  return mapping;
1532 }
1533 
1534 /**
1535  * @brief Create a reordering for hexahedra.
1536  */
1538 {
1539  const int order = conf.m_order;
1540  const int n = order - 1;
1541  const int n2 = n * n;
1542  int i, j, k;
1543 
1544  vector<int> mapping;
1545 
1546  // Map taking Gmsh edges to Nektar++ edges.
1547  static int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10};
1548 
1549  // Push back vertices.
1550  mapping.resize(8);
1551  for (i = 0; i < 8; ++i)
1552  {
1553  mapping[i] = i;
1554  }
1555 
1556  if (order == 1)
1557  {
1558  return mapping;
1559  }
1560 
1561  // Curvilinear edges
1562  mapping.resize(8 + 12 * n);
1563 
1564  // Reorder edges.
1565  int cnt = 8, offset;
1566  for (i = 0; i < 12; ++i)
1567  {
1568  int edge = gmshToNekEdge[i];
1569  offset = 8 + n * abs(edge);
1570 
1571  if (edge < 0)
1572  {
1573  for (int j = 0; j < n; ++j)
1574  {
1575  mapping[offset + n - j - 1] = cnt++;
1576  }
1577  }
1578  else
1579  {
1580  for (int j = 0; j < n; ++j)
1581  {
1582  mapping[offset + j] = cnt++;
1583  }
1584  }
1585  }
1586 
1587  if (conf.m_faceNodes == false || n2 == 0)
1588  {
1589  return mapping;
1590  }
1591 
1592  // Curvilinear face nodes.
1593  mapping.resize(8 + 12 * n + 6 * n2);
1594 
1595  // Map which takes Gmsh -> Nektar++ faces in the local element.
1596  static int gmsh2NekFace[6] = {0, 1, 4, 2, 3, 5};
1597 
1598  // Map which defines orientation between Gmsh and Nektar++ faces.
1599  StdRegions::Orientation faceOrient[6] = {
1605  StdRegions::eDir1FwdDir1_Dir2FwdDir2};
1606 
1607  for (i = 0; i < 6; ++i)
1608  {
1609  int face = gmsh2NekFace[i];
1610  int offset2 = 8 + 12 * n + i * n2;
1611  offset = 8 + 12 * n + face * n2;
1612 
1613  // Create a list of interior face nodes for this face only.
1614  vector<int> faceNodes(n2);
1615  for (j = 0; j < n2; ++j)
1616  {
1617  faceNodes[j] = offset2 + j;
1618  }
1619 
1620  // Now get the reordering of this face, which puts Gmsh
1621  // recursive ordering into Nektar++ row-by-row order.
1622  vector<int> tmp = quadTensorNodeOrdering(faceNodes, n);
1623 
1624  // Finally reorient the face according to the geometry
1625  // differences.
1626  if (faceOrient[i] == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
1627  {
1628  // Orientation is the same, just copy.
1629  for (j = 0; j < n2; ++j)
1630  {
1631  mapping[offset + j] = tmp[j];
1632  }
1633  }
1634  else if (faceOrient[i] == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
1635  {
1636  // Tranposed faces
1637  for (j = 0; j < n; ++j)
1638  {
1639  for (k = 0; k < n; ++k)
1640  {
1641  mapping[offset + j * n + k] = tmp[k * n + j];
1642  }
1643  }
1644  }
1645  else if (faceOrient[i] == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
1646  {
1647  for (j = 0; j < n; ++j)
1648  {
1649  for (k = 0; k < n; ++k)
1650  {
1651  mapping[offset + j * n + k] = tmp[j * n + (n - k - 1)];
1652  }
1653  }
1654  }
1655  }
1656 
1657  if (conf.m_volumeNodes == false)
1658  {
1659  return mapping;
1660  }
1661 
1662  const int totPoints = (order + 1) * (order + 1) * (order + 1);
1663  vector<int> interior;
1664  for (i = 8 + 12 * n + 6 * n2; i < totPoints; ++i)
1665  {
1666  interior.push_back(i);
1667  }
1668 
1669  interior = hexTensorNodeOrdering(interior, order - 1);
1670  mapping.insert(mapping.end(), interior.begin(), interior.end());
1671 
1672  return mapping;
1673 }
1674 
1675 /*
1676  * @brief Populate the element map #elmMap.
1677  *
1678  * This function primarily populates the element mapping #elmMap,
1679  * which takes a msh ID used by Gmsh and translates to element type,
1680  * element order and whether the element is incomplete (i.e. whether
1681  * it contains solely boundary nodes, or just face nodes). Note that
1682  * some of these elements, such as prisms of order >= 3, cannot yet be
1683  * generated by Gmsh.
1684  */
1685 std::map<unsigned int, ElmtConfig> InputGmsh::GenElmMap()
1686 {
1687  using namespace LibUtilities;
1688  std::map<unsigned int, ElmtConfig> tmp;
1689 
1690  // Elmt type, order, face, volume
1691  tmp[ 1] = ElmtConfig(eSegment, 1, false, false);
1692  tmp[ 2] = ElmtConfig(eTriangle, 1, false, false);
1693  tmp[ 3] = ElmtConfig(eQuadrilateral, 1, false, false);
1694  tmp[ 4] = ElmtConfig(eTetrahedron, 1, false, false);
1695  tmp[ 5] = ElmtConfig(eHexahedron, 1, false, false);
1696  tmp[ 6] = ElmtConfig(ePrism, 1, false, false);
1697  tmp[ 7] = ElmtConfig(ePyramid, 1, false, false);
1698  tmp[ 8] = ElmtConfig(eSegment, 2, true, false);
1699  tmp[ 9] = ElmtConfig(eTriangle, 2, false, false);
1700  tmp[ 10] = ElmtConfig(eQuadrilateral, 2, true, false);
1701  tmp[ 11] = ElmtConfig(eTetrahedron, 2, false, false);
1702  tmp[ 12] = ElmtConfig(eHexahedron, 2, true, true);
1703  tmp[ 13] = ElmtConfig(ePrism, 2, true, false);
1704  tmp[ 14] = ElmtConfig(ePyramid, 2, true, false);
1705  tmp[ 15] = ElmtConfig(ePoint, 1, false, false);
1706  tmp[ 16] = ElmtConfig(eQuadrilateral, 2, false, false);
1707  tmp[ 17] = ElmtConfig(eHexahedron, 2, false, false);
1708  tmp[ 18] = ElmtConfig(ePrism, 2, false, false);
1709  tmp[ 19] = ElmtConfig(ePyramid, 2, false, false);
1710  tmp[ 20] = ElmtConfig(eTriangle, 3, false, false);
1711  tmp[ 21] = ElmtConfig(eTriangle, 3, true, false);
1712  tmp[ 22] = ElmtConfig(eTriangle, 4, false, false);
1713  tmp[ 23] = ElmtConfig(eTriangle, 4, true, false);
1714  tmp[ 24] = ElmtConfig(eTriangle, 5, false, false);
1715  tmp[ 25] = ElmtConfig(eTriangle, 5, true, false);
1716  tmp[ 26] = ElmtConfig(eSegment, 3, true, false);
1717  tmp[ 27] = ElmtConfig(eSegment, 4, true, false);
1718  tmp[ 28] = ElmtConfig(eSegment, 5, true, false);
1719  tmp[ 29] = ElmtConfig(eTetrahedron, 3, true, false);
1720  tmp[ 30] = ElmtConfig(eTetrahedron, 4, true, true);
1721  tmp[ 31] = ElmtConfig(eTetrahedron, 5, true, true);
1722  tmp[ 32] = ElmtConfig(eTetrahedron, 4, true, false);
1723  tmp[ 33] = ElmtConfig(eTetrahedron, 5, true, false);
1724  tmp[ 36] = ElmtConfig(eQuadrilateral, 3, true, false);
1725  tmp[ 37] = ElmtConfig(eQuadrilateral, 4, true, false);
1726  tmp[ 38] = ElmtConfig(eQuadrilateral, 5, true, false);
1727  tmp[ 39] = ElmtConfig(eQuadrilateral, 3, false, false);
1728  tmp[ 40] = ElmtConfig(eQuadrilateral, 4, false, false);
1729  tmp[ 41] = ElmtConfig(eQuadrilateral, 5, false, false);
1730  tmp[ 42] = ElmtConfig(eTriangle, 6, true, false);
1731  tmp[ 43] = ElmtConfig(eTriangle, 7, true, false);
1732  tmp[ 44] = ElmtConfig(eTriangle, 8, true, false);
1733  tmp[ 45] = ElmtConfig(eTriangle, 9, true, false);
1734  tmp[ 46] = ElmtConfig(eTriangle, 10, true, false);
1735  tmp[ 47] = ElmtConfig(eQuadrilateral, 6, true, false);
1736  tmp[ 48] = ElmtConfig(eQuadrilateral, 7, true, false);
1737  tmp[ 49] = ElmtConfig(eQuadrilateral, 8, true, false);
1738  tmp[ 50] = ElmtConfig(eQuadrilateral, 9, true, false);
1739  tmp[ 51] = ElmtConfig(eQuadrilateral, 10, true, false);
1740  tmp[ 52] = ElmtConfig(eTriangle, 6, false, false);
1741  tmp[ 53] = ElmtConfig(eTriangle, 7, false, false);
1742  tmp[ 54] = ElmtConfig(eTriangle, 8, false, false);
1743  tmp[ 55] = ElmtConfig(eTriangle, 9, false, false);
1744  tmp[ 56] = ElmtConfig(eTriangle, 10, false, false);
1745  tmp[ 57] = ElmtConfig(eQuadrilateral, 6, false, false);
1746  tmp[ 58] = ElmtConfig(eQuadrilateral, 7, false, false);
1747  tmp[ 59] = ElmtConfig(eQuadrilateral, 8, false, false);
1748  tmp[ 60] = ElmtConfig(eQuadrilateral, 9, false, false);
1749  tmp[ 61] = ElmtConfig(eQuadrilateral, 10, false, false);
1750  tmp[ 62] = ElmtConfig(eSegment, 6, true, false);
1751  tmp[ 63] = ElmtConfig(eSegment, 7, true, false);
1752  tmp[ 64] = ElmtConfig(eSegment, 8, true, false);
1753  tmp[ 65] = ElmtConfig(eSegment, 9, true, false);
1754  tmp[ 66] = ElmtConfig(eSegment, 10, true, false);
1755  tmp[ 71] = ElmtConfig(eTetrahedron, 6, true, true);
1756  tmp[ 72] = ElmtConfig(eTetrahedron, 7, true, true);
1757  tmp[ 73] = ElmtConfig(eTetrahedron, 8, true, true);
1758  tmp[ 74] = ElmtConfig(eTetrahedron, 9, true, true);
1759  tmp[ 75] = ElmtConfig(eTetrahedron, 10, true, true);
1760  tmp[ 79] = ElmtConfig(eTetrahedron, 6, true, false);
1761  tmp[ 80] = ElmtConfig(eTetrahedron, 7, true, false);
1762  tmp[ 81] = ElmtConfig(eTetrahedron, 8, true, false);
1763  tmp[ 82] = ElmtConfig(eTetrahedron, 9, true, false);
1764  tmp[ 83] = ElmtConfig(eTetrahedron, 10, true, false);
1765  tmp[ 90] = ElmtConfig(ePrism, 3, true, true);
1766  tmp[ 91] = ElmtConfig(ePrism, 4, true, true);
1767  tmp[ 92] = ElmtConfig(eHexahedron, 3, true, true);
1768  tmp[ 93] = ElmtConfig(eHexahedron, 4, true, true);
1769  tmp[ 94] = ElmtConfig(eHexahedron, 5, true, true);
1770  tmp[ 95] = ElmtConfig(eHexahedron, 6, true, true);
1771  tmp[ 96] = ElmtConfig(eHexahedron, 7, true, true);
1772  tmp[ 97] = ElmtConfig(eHexahedron, 8, true, true);
1773  tmp[ 98] = ElmtConfig(eHexahedron, 9, true, true);
1774  tmp[ 99] = ElmtConfig(eHexahedron, 3, true, false);
1775  tmp[100] = ElmtConfig(eHexahedron, 4, true, false);
1776  tmp[101] = ElmtConfig(eHexahedron, 5, true, false);
1777  tmp[102] = ElmtConfig(eHexahedron, 6, true, false);
1778  tmp[103] = ElmtConfig(eHexahedron, 7, true, false);
1779  tmp[104] = ElmtConfig(eHexahedron, 8, true, false);
1780  tmp[105] = ElmtConfig(eHexahedron, 9, true, false);
1781  tmp[106] = ElmtConfig(ePrism, 5, true, true);
1782  tmp[107] = ElmtConfig(ePrism, 6, true, true);
1783  tmp[108] = ElmtConfig(ePrism, 7, true, true);
1784  tmp[109] = ElmtConfig(ePrism, 8, true, true);
1785  tmp[110] = ElmtConfig(ePrism, 9, true, true);
1786  tmp[111] = ElmtConfig(ePrism, 3, true, false);
1787  tmp[112] = ElmtConfig(ePrism, 4, true, false);
1788  tmp[113] = ElmtConfig(ePrism, 5, true, false);
1789  tmp[114] = ElmtConfig(ePrism, 6, true, false);
1790  tmp[115] = ElmtConfig(ePrism, 7, true, false);
1791  tmp[116] = ElmtConfig(ePrism, 8, true, false);
1792  tmp[117] = ElmtConfig(ePrism, 9, true, false);
1793  tmp[118] = ElmtConfig(ePyramid, 3, true, true);
1794  tmp[119] = ElmtConfig(ePyramid, 4, true, true);
1795  tmp[120] = ElmtConfig(ePyramid, 5, true, true);
1796  tmp[121] = ElmtConfig(ePyramid, 6, true, true);
1797  tmp[122] = ElmtConfig(ePyramid, 7, true, true);
1798  tmp[123] = ElmtConfig(ePyramid, 8, true, true);
1799  tmp[124] = ElmtConfig(ePyramid, 9, true, true);
1800  tmp[125] = ElmtConfig(ePyramid, 3, true, false);
1801  tmp[126] = ElmtConfig(ePyramid, 4, true, false);
1802  tmp[127] = ElmtConfig(ePyramid, 5, true, false);
1803  tmp[128] = ElmtConfig(ePyramid, 6, true, false);
1804  tmp[129] = ElmtConfig(ePyramid, 7, true, false);
1805  tmp[130] = ElmtConfig(ePyramid, 7, true, false);
1806  tmp[131] = ElmtConfig(ePyramid, 8, true, false);
1807 
1808  return tmp;
1809 }
1810 
1811 }
1812 }
bool m_faceNodes
Denotes whether the element contains face nodes. For 2D elements, if this is true then the element co...
Definition: ElementConfig.h:80
std::vector< int > prismTensorNodeOrdering(const std::vector< int > &nodes, int n)
Reorder Gmsh nodes so that they appear in a "tensor product" format suitable for the interior of a Ne...
Definition: InputGmsh.cpp:464
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
std::vector< int > triTensorNodeOrdering(const std::vector< int > &nodes, int n)
Reorder Gmsh nodes so that they appear in a "tensor product" format suitable for the interior of a Ne...
Definition: InputGmsh.cpp:170
Basic information about an element.
Definition: ElementConfig.h:50
io::filtering_istream m_mshFile
Input stream.
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: HOAlignment.h:135
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 Gmsh nodes so that they appear in a tensor product format suitable for the interior of a Nekt...
Definition: InputGmsh.cpp:93
std::vector< T > surfVerts
The quadrilateral surface vertices – templated so that this can either be nodes or IDs...
Definition: HOAlignment.h:231
STL namespace.
std::vector< T > surfVerts
The triangle surface vertices – templated so that this can either be nodes or IDs.
Definition: HOAlignment.h:65
pair< ModuleType, string > ModuleKey
static std::vector< int > PrismReordering(NekMeshUtils::ElmtConfig conf)
Create a reordering for prisms.
Definition: InputGmsh.cpp:1318
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
static std::vector< int > TetReordering(NekMeshUtils::ElmtConfig conf)
Create a reordering for tetrahedra.
Definition: InputGmsh.cpp:1166
std::vector< int > hexTensorNodeOrdering(const std::vector< int > &nodes, int n)
Reorder Gmsh nodes so that they appear in a "tensor product" format suitable for the interior of a Ne...
Definition: InputGmsh.cpp:540
const char *const ShapeTypeMap[]
Definition: ShapeType.hpp:66
unsigned int m_order
Order of the element.
Definition: ElementConfig.h:87
static std::vector< int > HexReordering(NekMeshUtils::ElmtConfig conf)
Create a reordering for hexahedra.
Definition: InputGmsh.cpp:1537
NEKMESHUTILS_EXPORT void OpenStream()
Open a file for input.
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a line.
Definition: Line.cpp:187
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a prism.
Definition: Prism.cpp:234
static std::vector< int > CreateReordering(unsigned int InputGmshEntity)
Create a reordering map for a given element.
Definition: InputGmsh.cpp:1011
bool m_volumeNodes
Denotes whether the element contains volume (i.e. interior) nodes. These are not supported by either ...
Definition: ElementConfig.h:85
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
static std::vector< int > LineReordering(NekMeshUtils::ElmtConfig conf)
Definition: InputGmsh.cpp:1053
bool operator()(Mode const &a, Mode const &b) const
Definition: InputGmsh.cpp:224
static std::vector< int > QuadReordering(NekMeshUtils::ElmtConfig conf)
Create a reordering for quadrilaterals.
Definition: InputGmsh.cpp:1118
void Align(std::vector< int > vertId)
Align this surface to a given vertex ID.
Definition: HOAlignment.h:277
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a triangle.
A lightweight struct for dealing with high-order triangle alignment.
Definition: HOAlignment.h:50
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a hexahedron.
Definition: Hexahedron.cpp:271
Abstract base class for input modules.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
boost::tuple< int, int, int > Mode
Definition: InputGmsh.cpp:221
A lightweight struct for dealing with high-order quadrilateral alignment.
Definition: HOAlignment.h:215
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a tetrahedron.
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
std::vector< int > tetTensorNodeOrdering(const std::vector< int > &nodes, int n)
Reorder Gmsh nodes so that they appear in a "tensor product" format suitable for the interior of a Ne...
Definition: InputGmsh.cpp:285
int GetNnodes(unsigned int InputGmshEntity)
Definition: InputGmsh.cpp:954
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a pyramid.
Definition: Pyramid.cpp:175
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
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:72
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
static std::vector< int > TriReordering(NekMeshUtils::ElmtConfig conf)
Create a reordering for triangles.
Definition: InputGmsh.cpp:1069
std::pair< ModuleType, std::string > ModuleKey
static std::map< unsigned int, NekMeshUtils::ElmtConfig > GenElmMap()
Definition: InputGmsh.cpp:1685
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
static std::map< unsigned int, NekMeshUtils::ElmtConfig > elmMap
Definition: InputGmsh.h:68
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215