MFEM  v4.4.0 Finite element discretization library
gmsh.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "gmsh.hpp"
13 #include "vtk.hpp"
14
15 namespace mfem
16 {
17
18 int BarycentricToGmshTet(int *b, int ref)
19 {
20  int i = b[0];
21  int j = b[1];
22  int k = b[2];
23  int l = b[3];
24  bool ibdr = (i == 0);
25  bool jbdr = (j == 0);
26  bool kbdr = (k == 0);
27  bool lbdr = (l == 0);
28  if (ibdr && jbdr && kbdr)
29  {
30  return 0;
31  }
32  else if (jbdr && kbdr && lbdr)
33  {
34  return 1;
35  }
36  else if (ibdr && kbdr && lbdr)
37  {
38  return 2;
39  }
40  else if (ibdr && jbdr && lbdr)
41  {
42  return 3;
43  }
44  int offset = 4;
45  if (jbdr && kbdr) // Edge DOF on j == 0 and k == 0
46  {
47  return offset + i - 1;
48  }
49  else if (kbdr && lbdr) // Edge DOF on k == 0 and l == 0
50  {
51  return offset + ref - 1 + j - 1;
52  }
53  else if (ibdr && kbdr) // Edge DOF on i == 0 and k == 0
54  {
55  return offset + 2 * (ref - 1) + ref - j - 1;
56  }
57  else if (ibdr && jbdr) // Edge DOF on i == 0 and j == 0
58  {
59  return offset + 3 * (ref - 1) + ref - k - 1;
60  }
61  else if (ibdr && lbdr) // Edge DOF on i == 0 and l == 0
62  {
63  return offset + 4 * (ref - 1) + ref - k - 1;
64  }
65  else if (jbdr && lbdr) // Edge DOF on j == 0 and l == 0
66  {
67  return offset + 5 * (ref - 1) + ref - k - 1;
68  }
69
70  // Recursive numbering for the faces
71  offset += 6 * (ref - 1);
72  if (kbdr)
73  {
74  int b_out[3];
75  b_out[0] = j-1;
76  b_out[1] = i-1;
77  b_out[2] = ref - i - j - 1;
78  return offset + BarycentricToVTKTriangle(b_out, ref-3);
79  }
80  else if (jbdr)
81  {
82  int b_out[3];
83  b_out[0] = i-1;
84  b_out[1] = k-1;
85  b_out[2] = ref - i - k - 1;
86  offset += (ref - 1) * (ref - 2) / 2;
87  return offset + BarycentricToVTKTriangle(b_out, ref-3);
88  }
89  else if (ibdr)
90  {
91  int b_out[3];
92  b_out[0] = k-1;
93  b_out[1] = j-1;
94  b_out[2] = ref - j - k - 1;
95  offset += (ref - 1) * (ref - 2);
96  return offset + BarycentricToVTKTriangle(b_out, ref-3);
97  }
98  else if (lbdr)
99  {
100  int b_out[3];
101  b_out[0] = ref-j-k-1;
102  b_out[1] = j-1;
103  b_out[2] = k-1;
104  offset += 3 * (ref - 1) * (ref - 2) / 2;
105  return offset + BarycentricToVTKTriangle(b_out, ref-3);
106  }
107
108  // Recursive numbering for interior
109  {
110  int b_out[4];
111  b_out[0] = i-1;
112  b_out[1] = j-1;
113  b_out[2] = k-1;
114  b_out[3] = ref - i - j - k - 1;
115  offset += 2 * (ref - 1) * (ref - 2);
116  return offset + BarycentricToGmshTet(b_out, ref-4);
117  }
118 }
119
120 int CartesianToGmshQuad(int idx_in[], int ref)
121 {
122  int i = idx_in[0];
123  int j = idx_in[1];
124  // Do we lie on any of the edges
125  bool ibdr = (i == 0 || i == ref);
126  bool jbdr = (j == 0 || j == ref);
127  if (ibdr && jbdr) // Vertex DOF
128  {
129  return (i ? (j ? 2 : 1) : (j ? 3 : 0));
130  }
131  int offset = 4;
132  if (jbdr) // Edge DOF on j==0 or j==ref
133  {
134  return offset + (j ? 3*ref - 3 - i : i - 1);
135  }
136  else if (ibdr) // Edge DOF on i==0 or i==ref
137  {
138  return offset + (i ? ref - 1 + j - 1 : 4*ref - 4 - j);
139  }
140  else // Recursive numbering for interior
141  {
142  int idx_out[2];
143  idx_out[0] = i-1;
144  idx_out[1] = j-1;
145  offset += 4 * (ref - 1);
146  return offset + CartesianToGmshQuad(idx_out, ref-2);
147  }
148 }
149
150 int CartesianToGmshHex(int idx_in[], int ref)
151 {
152  int i = idx_in[0];
153  int j = idx_in[1];
154  int k = idx_in[2];
155  // Do we lie on any of the edges
156  bool ibdr = (i == 0 || i == ref);
157  bool jbdr = (j == 0 || j == ref);
158  bool kbdr = (k == 0 || k == ref);
159  if (ibdr && jbdr && kbdr) // Vertex DOF
160  {
161  return (i ? (j ? (k ? 6 : 2) : (k ? 5 : 1)) :
162  (j ? (k ? 7 : 3) : (k ? 4 : 0)));
163  }
164  int offset = 8;
165  if (jbdr && kbdr) // Edge DOF on x-directed edge
166  {
167  return offset + (j ? (k ? 12*ref-12-i: 6*ref-6-i) :
168  (k ? 8*ref-9+i: i-1));
169  }
170  else if (ibdr && kbdr) // Edge DOF on y-directed edge
171  {
172  return offset + (k ? (i ? 10*ref-11+j: 9*ref-10+j) :
173  (i ? 3*ref-4+j: ref-2+j));
174  }
175  else if (ibdr && jbdr) // Edge DOF on z-directed edge
176  {
177  return offset + (i ? (j ? 6*ref-7+k: 4*ref-5+k) :
178  (j ? 7*ref-8+k: 2*ref-3+k));
179  }
180  else if (ibdr) // Face DOF on x-directed face
181  {
182  int idx_out[2];
183  idx_out[0] = i ? j-1 : k-1;
184  idx_out[1] = i ? k-1 : j-1;
185  offset += (12 + (i ? 3 : 2) * (ref - 1)) * (ref - 1);
186  return offset + CartesianToGmshQuad(idx_out, ref-2);
187  }
188  else if (jbdr) // Face DOF on y-directed face
189  {
190  int idx_out[2];
191  idx_out[0] = j ? ref-i-1 : i-1;
192  idx_out[1] = j ? k-1 : k-1;
193  offset += (12 + (j ? 4 : 1) * (ref - 1)) * (ref - 1);
194  return offset + CartesianToGmshQuad(idx_out, ref-2);
195  }
196  else if (kbdr) // Face DOF on z-directed face
197  {
198  int idx_out[2];
199  idx_out[0] = k ? i-1 : j-1;
200  idx_out[1] = k ? j-1 : i-1;
201  offset += (12 + (k ? 5 : 0) * (ref - 1)) * (ref - 1);
202  return offset + CartesianToGmshQuad(idx_out, ref-2);
203  }
204  else // Recursive numbering for interior
205  {
206  int idx_out[3];
207  idx_out[0] = i-1;
208  idx_out[1] = j-1;
209  idx_out[2] = k-1;
210
211  offset += (12 + 6 * (ref - 1)) * (ref - 1);
212  return offset + CartesianToGmshHex(idx_out, ref-2);
213  }
214 }
215
216 int WedgeToGmshPri(int idx_in[], int ref)
217 {
218  int i = idx_in[0];
219  int j = idx_in[1];
220  int k = idx_in[2];
221  int l = ref - i -j;
222  bool ibdr = (i == 0);
223  bool jbdr = (j == 0);
224  bool kbdr = (k == 0 || k == ref);
225  bool lbdr = (l == 0);
226  if (ibdr && jbdr && kbdr)
227  {
228  return k ? 3 : 0;
229  }
230  else if (jbdr && lbdr && kbdr)
231  {
232  return k ? 4 : 1;
233  }
234  else if (ibdr && lbdr && kbdr)
235  {
236  return k ? 5 : 2;
237  }
238  int offset = 6;
239  if (jbdr && kbdr)
240  {
241  return offset + (k ? 6 * (ref - 1) + i - 1: i - 1);
242  }
243  else if (ibdr && kbdr)
244  {
245  return offset + (k ? 7 * (ref -1) + j-1 : ref - 1 + j - 1);
246  }
247  else if (ibdr && jbdr)
248  {
249  return offset + 2 * (ref - 1) + k - 1;
250  }
251  else if (lbdr && kbdr)
252  {
253  return offset + (k ? 8 * (ref -1) + j - 1 : 3 * (ref - 1) + j - 1);
254  }
255  else if (jbdr && lbdr)
256  {
257  return offset + 4 * (ref - 1) + k - 1;
258  }
259  else if (ibdr && lbdr)
260  {
261  return offset + 5 * (ref - 1) + k - 1;
262  }
263  offset += 9 * (ref-1);
264  if (kbdr) // Triangular faces at k=0 and k=ref
265  {
266  int b_out[3];
267  b_out[0] = k ? i-1 : j-1;
268  b_out[1] = k ? j-1 : i-1;
269  b_out[2] = ref - i - j - 1;
270  offset += k ? (ref-1)*(ref-2) / 2: 0;
271  return offset + BarycentricToVTKTriangle(b_out, ref-3);
272  }
273  offset += (ref-1)*(ref-2);
274  if (jbdr) // Quadrilateral face at j=0
275  {
276  int idx_out[2];
277  idx_out[0] = i-1;
278  idx_out[1] = k-1;
279  return offset + CartesianToGmshQuad(idx_out, ref-2);
280  }
281  else if (ibdr) // Quadrilateral face at i=0
282  {
283  int idx_out[2];
284  idx_out[0] = k-1;
285  idx_out[1] = j-1;
286  offset += (ref-1)*(ref-1);
287  return offset + CartesianToGmshQuad(idx_out, ref-2);
288  }
289  else if (lbdr) // Quadrilateral face at l=ref-i-j=0
290  {
291  int idx_out[2];
292  idx_out[0] = j-1;
293  idx_out[1] = k-1;
294  offset += 2*(ref-1)*(ref-1);
295  return offset + CartesianToGmshQuad(idx_out, ref-2);
296  }
297  offset += 3*(ref-1)*(ref-1);
298  // The Gmsh Prism interiors are a tensor product of segments of order ref-2
299  // and triangles of order ref-3
300  {
301  int b_out[3];
302  b_out[0] = i-1;
303  b_out[1] = j-1;
304  b_out[2] = ref - i - j - 1;
305  int ot = BarycentricToVTKTriangle(b_out, ref-3);
306  int os = (k==1) ? 0 : (k == ref-1 ? 1 : k);
307  return offset + (ref-1) * ot + os;
308  }
309 }
310
311 int CartesianToGmshPyramid(int idx_in[], int ref)
312 {
313  int i = idx_in[0];
314  int j = idx_in[1];
315  int k = idx_in[2];
316  // Do we lie on any of the edges
317  bool ibdr = (i == 0 || i == ref-k);
318  bool jbdr = (j == 0 || j == ref-k);
319  bool kbdr = (k == 0);
320  if (ibdr && jbdr && kbdr)
321  {
322  return i ? (j ? 2 : 1): (j ? 3 : 0);
323  }
324  else if (k == ref)
325  {
326  return 4;
327  }
328  int offset = 5;
329  if (jbdr && kbdr)
330  {
331  return offset + (j ? (6 * ref - 6 - i) : (i - 1));
332  }
333  else if (ibdr && kbdr)
334  {
335  return offset + (i ? (3 * ref - 4 + j) : (ref - 2 + j));
336  }
337  else if (ibdr && jbdr)
338  {
339  return offset + (i ? (j ? 6 : 4) : (j ? 7 : 2 )) * (ref-1) + k - 1;
340  }
341  offset += 8*(ref-1);
342  if (jbdr)
343  {
344  int b_out[3];
345  b_out[0] = j ? ref - i - k - 1 : i - 1;
346  b_out[1] = k - 1;
347  b_out[2] = (j ? i - 1 : ref - i - k - 1);
348  offset += (j ? 3 : 0) * (ref - 1) * (ref - 2) / 2;
349  return offset + BarycentricToVTKTriangle(b_out, ref-3);
350  }
351  else if (ibdr)
352  {
353  int b_out[3];
354  b_out[0] = i ? j - 1: ref - j - k - 1;
355  b_out[1] = k - 1;
356  b_out[2] = (i ? ref - j - k - 1: j - 1);
357  offset += (i ? 2 : 1) * (ref - 1) * (ref - 2) / 2;
358  return offset + BarycentricToVTKTriangle(b_out, ref-3);
359  }
360  else if (kbdr)
361  {
362  int idx_out[2];
363  idx_out[0] = k ? i-1 : j-1;
364  idx_out[1] = k ? j-1 : i-1;
365  offset += 2 * (ref - 1) * (ref - 2);
366  return offset + CartesianToGmshQuad(idx_out, ref-2);
367  }
368  offset += (2 * (ref - 2) + (ref - 1)) * (ref - 1) ;
369  {
370  int idx_out[3];
371  idx_out[0] = i-1;
372  idx_out[1] = j-1;
373  idx_out[2] = k-1;
374  return offset + CartesianToGmshPyramid(idx_out, ref-3);
375  }
376 }
377
378 void GmshHOSegmentMapping(int order, int *map)
379 {
380  map[0] = 0;
381  map[order] = 1;
382  for (int i=1; i<order; i++)
383  {
384  map[i] = i + 1;
385  }
386 }
387
388 void GmshHOTriangleMapping(int order, int *map)
389 {
390  int b[3];
391  int o = 0;
392  for (b[1]=0; b[1]<=order; ++b[1])
393  {
394  for (b[0]=0; b[0]<=order-b[1]; ++b[0])
395  {
396  b[2] = order - b[0] - b[1];
397  map[o] = BarycentricToVTKTriangle(b, order);
398  o++;
399  }
400  }
401 }
402
403 void GmshHOQuadrilateralMapping(int order, int *map)
404 {
405  int b[2];
406  int o = 0;
407  for (b[1]=0; b[1]<=order; b[1]++)
408  {
409  for (b[0]=0; b[0]<=order; b[0]++)
410  {
412  o++;
413  }
414  }
415 }
416
417 void GmshHOTetrahedronMapping(int order, int *map)
418 {
419  int b[4];
420  int o = 0;
421  for (b[2]=0; b[2]<=order; ++b[2])
422  {
423
424  for (b[1]=0; b[1]<=order-b[2]; ++b[1])
425  {
426  for (b[0]=0; b[0]<=order-b[1]-b[2]; ++b[0])
427  {
428  b[3] = order - b[0] - b[1] - b[2];
429  map[o] = BarycentricToGmshTet(b, order);
430  o++;
431  }
432  }
433  }
434 }
435
436 void GmshHOHexahedronMapping(int order, int *map)
437 {
438  int b[3];
439  int o = 0;
440  for (b[2]=0; b[2]<=order; b[2]++)
441  {
442  for (b[1]=0; b[1]<=order; b[1]++)
443  {
444  for (b[0]=0; b[0]<=order; b[0]++)
445  {
446  map[o] = CartesianToGmshHex(b, order);
447  o++;
448  }
449  }
450  }
451 }
452
453 void GmshHOWedgeMapping(int order, int *map)
454 {
455  int b[3];
456  int o = 0;
457  for (b[2]=0; b[2]<=order; b[2]++)
458  {
459  for (b[1]=0; b[1]<=order; b[1]++)
460  {
461  for (b[0]=0; b[0]<=order - b[1]; b[0]++)
462  {
463  map[o] = WedgeToGmshPri(b, order);
464  o++;
465  }
466  }
467  }
468 }
469
470 void GmshHOPyramidMapping(int order, int *map)
471 {
472  int b[3];
473  int o = 0;
474  for (b[2]=0; b[2]<=order; b[2]++)
475  {
476  for (b[1]=0; b[1]<=order - b[2]; b[1]++)
477  {
478  for (b[0]=0; b[0]<=order - b[2]; b[0]++)
479  {
480  map[o] = CartesianToGmshPyramid(b, order);
481  o++;
482  }
483  }
484  }
485 }
486
487 } // namespace mfem
void GmshHOSegmentMapping(int order, int *map)
Generate Gmsh vertex mapping for a Segment.
Definition: gmsh.cpp:378
int CartesianToGmshHex(int idx_in[], int ref)
Definition: gmsh.cpp:150
int BarycentricToGmshTet(int *b, int ref)
Definition: gmsh.cpp:18
int BarycentricToVTKTriangle(int *b, int ref)
Return the VTK node index of the barycentric point b in a triangle with refinement level ref...
Definition: vtk.cpp:153
double b
Definition: lissajous.cpp:42
Generate Gmsh vertex mapping for a Quadrilateral.
Definition: gmsh.cpp:403
void GmshHOTriangleMapping(int order, int *map)
Generate Gmsh vertex mapping for a Triangle.
Definition: gmsh.cpp:388
void GmshHOTetrahedronMapping(int order, int *map)
Generate Gmsh vertex mapping for a Tetrahedron.
Definition: gmsh.cpp:417
int WedgeToGmshPri(int idx_in[], int ref)
Definition: gmsh.cpp:216