1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributor(s):
7 // Thomas Toulorge
8
9 #include <iostream>
10 #include <vector>
11 #include <set>
12 #include <algorithm>
13 #include "SVector3.h"
14 #include "GmshDefines.h"
15 #include "fullMatrix.h"
16 #include "ElementType.h"
17 #include "BasisFactory.h"
18 #include "nodalBasis.h"
19 #include "CGNSCommon.h"
20 #include "CGNSConventions.h"
21
22 #if defined(HAVE_LIBCGNS)
23
24 namespace {
25
26 // ------------- Conversion of node ordering between CGNS and Gmsh
27 // -------------
28
addEdgePointsCGNS(const SVector3 p0,const SVector3 p1,int order,std::vector<SVector3> & points)29 void addEdgePointsCGNS(const SVector3 p0, const SVector3 p1, int order,
30 std::vector<SVector3> &points)
31 {
32 if(order < 2) return;
33
34 double ds = 1. / order;
35 for(int i = 1; i < order; i++) {
36 double f = ds * i;
37 points.push_back(p0 * (1. - f) + p1 * f);
38 }
39 }
40
41 std::vector<SVector3> generatePointsTriCGNS(int order, bool complete);
42
addTriPointsCGNS(const SVector3 p0,const SVector3 p1,const SVector3 p2,int order,std::vector<SVector3> & points)43 void addTriPointsCGNS(const SVector3 p0, const SVector3 p1, const SVector3 p2,
44 int order, std::vector<SVector3> &points)
45 {
46 if(order < 3) return;
47
48 std::vector<SVector3> triPoints = generatePointsTriCGNS(order - 3, true);
49
50 double scale = double(order - 3) / double(order);
51 SVector3 offset(1. / order, 1. / order, 0);
52
53 for(size_t i = 0; i < triPoints.size(); i++) {
54 SVector3 ip = triPoints[i];
55 double u = ip[0] * scale + 1. / order;
56 double v = ip[1] * scale + 1. / order;
57 SVector3 pt = (1. - u - v) * p0 + u * p1 + v * p2;
58 points.push_back(pt);
59 }
60 }
61
addQuaPointsCGNS(int order,std::vector<SVector3> & points)62 void addQuaPointsCGNS(int order, std::vector<SVector3> &points)
63 {
64 if(order > 2) {
65 double scale = double(order - 2) / double(order);
66
67 SVector3 corner[4] = {SVector3(-1, -1, 0), SVector3(1, -1, 0),
68 SVector3(1, 1, 0), SVector3(-1, 1, 0)};
69
70 for(int i = 0; i < 4; i++) {
71 SVector3 c1 = corner[i];
72 SVector3 c2 = corner[(i + 1) % 4];
73 double ds = 1. / (order - 2);
74 for(int i = 0; i < order - 2; i++)
75 points.push_back((c1 * (1. - i * ds) + c2 * (i * ds)) * scale);
76 }
77 }
78
79 if(order == 2 || order == 4) points.push_back(SVector3(0, 0, 0));
80 }
81
addQuaPointsCGNS(const SVector3 p0,const SVector3 p1,const SVector3 p2,const SVector3 p3,int order,std::vector<SVector3> & points)82 void addQuaPointsCGNS(const SVector3 p0, const SVector3 p1, const SVector3 p2,
83 const SVector3 p3, int order,
84 std::vector<SVector3> &points)
85 {
86 std::vector<SVector3> quaPoints;
87
88 addQuaPointsCGNS(order, quaPoints);
89
90 for(size_t i = 0; i < quaPoints.size(); i++) {
91 SVector3 ip = quaPoints[i];
92 double u = ip[0];
93 double v = ip[1];
94 SVector3 pt = ((1. - u) * (1. - v) * p0 + (1. + u) * (1. - v) * p1 +
95 (1. + u) * (1. + v) * p2 + (1. - u) * (1. + v) * p3) *
96 0.25;
97 points.push_back(pt);
98 }
99 }
100
101 // void print(std::vector<SVector3> &points, const char *title, int iStart,
102 // int iEnd = -1)
103 // {
104 // iEnd = iEnd == -1 ? points.size() : iEnd;
105
106 // std::cout << title << std::endl;
107 // for(int i = iStart; i < iEnd; i++) {
108 // std::cout << i << " :";
109 // for(int d = 0; d < 3; d++) std::cout << " " << points[i][d];
110 // std::cout << std::endl;
111 // }
112 // }
113
generatePointsEdgeCGNS(int order)114 std::vector<SVector3> generatePointsEdgeCGNS(int order)
115 {
116 std::vector<SVector3> pp;
117
118 if(order == 0) {
119 pp.push_back(SVector3(0., 0., 0.));
120 return pp;
121 }
122
123 // primary vertices
124 pp.push_back(SVector3(-1, 0, 0));
125 pp.push_back(SVector3(1, 0, 0));
126
127 // internal points
128 addEdgePointsCGNS(pp[0], pp[1], order, pp);
129
130 return pp;
131 }
132
generatePointsTriCGNS(int order,bool complete)133 std::vector<SVector3> generatePointsTriCGNS(int order, bool complete)
134 {
135 std::vector<SVector3> pp;
136
137 if(order == 0) {
138 pp.push_back(SVector3(1. / 3., 1. / 3., 0.));
139 return pp;
140 }
141
142 // primary vertices
143 pp.push_back(SVector3(0, 0, 0));
144 pp.push_back(SVector3(1, 0, 0));
145 pp.push_back(SVector3(0, 1, 0));
146
147 // internal points of edges
148 for(int i = 0; i < 3; i++) {
149 addEdgePointsCGNS(pp[i], pp[(i + 1) % 3], order, pp);
150 }
151
152 // internal points
153 if(complete && order > 2) {
154 addTriPointsCGNS(pp[0], pp[1], pp[2], order, pp);
155 }
156
157 return pp;
158 }
159
generatePointsQuaCGNS(int order,bool complete)160 std::vector<SVector3> generatePointsQuaCGNS(int order, bool complete)
161 {
162 std::vector<SVector3> pp;
163
164 if(order == 0) {
165 pp.push_back(SVector3(0., 0., 0.));
166 return pp;
167 }
168
169 // primary vertices
170 pp.push_back(SVector3(-1, -1, 0));
171 pp.push_back(SVector3(1, -1, 0));
172 pp.push_back(SVector3(1, 1, 0));
173 pp.push_back(SVector3(-1, 1, 0));
174
175 // internal points of edges
176 for(int i = 0; i < 4; i++) {
177 addEdgePointsCGNS(pp[i], pp[(i + 1) % 4], order, pp);
178 }
179
180 // internal points
181 if(complete && order > 1) {
182 addQuaPointsCGNS(pp[0], pp[1], pp[2], pp[3], order, pp);
183 }
184
185 return pp;
186 }
187
generatePointsTetCGNS(int order,bool complete)188 std::vector<SVector3> generatePointsTetCGNS(int order, bool complete)
189 {
190 std::vector<SVector3> pp;
191
192 if(order == 0) {
193 pp.push_back(SVector3(1. / 4., 1. / 4., 1. / 4.));
194 return pp;
195 }
196
197 // primary vertices
198 pp.push_back(SVector3(0, 0, 0));
199 pp.push_back(SVector3(1, 0, 0));
200 pp.push_back(SVector3(0, 1, 0));
201 pp.push_back(SVector3(0, 0, 1));
202
203 // internal points in edges of base triangle
204 for(int i = 0; i < 3; i++)
205 addEdgePointsCGNS(pp[i], pp[(i + 1) % 3], order, pp);
206
207 // internal points in upstanding edges
208 for(int i = 0; i < 3; i++) addEdgePointsCGNS(pp[i], pp[3], order, pp);
209
210 if(complete && order > 2) {
211 // internal points of base triangle
212 addTriPointsCGNS(pp[0], pp[1], pp[2], order, pp);
213
214 // internal points of upstanding triangles
215 for(int i = 0; i < 3; i++) {
216 addTriPointsCGNS(pp[i], pp[(i + 1) % 3], pp[3], order, pp);
217 }
218
219 // internal points as a tet of order p-3
220 if(order > 3) {
221 std::vector<SVector3> tetPp = generatePointsTetCGNS(order - 4, true);
222
223 double scale = (order - 4) / order;
224 SVector3 offset(1. / order, 1. / order, 1. / order);
225 for(size_t i = 0; i < tetPp.size(); i++) {
226 SVector3 volumePoint = tetPp[i];
227 volumePoint *= scale;
228 volumePoint += offset;
229 pp.push_back(volumePoint);
230 }
231 }
232 }
233
234 return pp;
235 }
236
generatePointsHexCGNS(int order,bool complete)237 std::vector<SVector3> generatePointsHexCGNS(int order, bool complete)
238 {
239 std::vector<SVector3> pp;
240
241 if(order == 0) {
242 pp.push_back(SVector3(0., 0., 0.));
243 return pp;
244 }
245
246 // principal vertices
247 pp.push_back(SVector3(-1, -1, -1));
248 pp.push_back(SVector3(1, -1, -1));
249 pp.push_back(SVector3(1, 1, -1));
250 pp.push_back(SVector3(-1, 1, -1));
251 pp.push_back(SVector3(-1, -1, 1));
252 pp.push_back(SVector3(1, -1, 1));
253 pp.push_back(SVector3(1, 1, 1));
254 pp.push_back(SVector3(-1, 1, 1));
255
256 // internal points of base quadrangle edges
257 for(int i = 0; i < 4; i++)
258 addEdgePointsCGNS(pp[i], pp[(i + 1) % 4], order, pp);
259
260 std::vector<SVector3> up[4];
261 // internal points of mounting edges
262 for(int i = 0; i < 4; i++) {
263 addEdgePointsCGNS(pp[i], pp[i + 4], order, up[i]);
264 pp.insert(pp.end(), up[i].begin(), up[i].end());
265 }
266
267 // internal points of top quadrangle edges
268 for(int i = 0; i < 4; i++)
269 addEdgePointsCGNS(pp[i + 4], pp[(i + 1) % 4 + 4], order, pp);
270
271 if(complete && order > 1) {
272 // internal points of base quadrangle
273 addQuaPointsCGNS(pp[0], pp[1], pp[2], pp[3], order, pp);
274
275 // internal points of upstanding faces
276 for(int i = 0; i < 4; i++) {
277 addQuaPointsCGNS(pp[i], pp[(i + 1) % 4], pp[(i + 1) % 4 + 4], pp[i + 4],
278 order, pp);
279 }
280
281 // internal points of top quadrangle
282 addQuaPointsCGNS(pp[4], pp[5], pp[6], pp[7], order, pp);
283
284 // internal volume points as a succession of internal planes
285 for(int i = 0; i <= order - 2; i++) {
286 addQuaPointsCGNS(up[0][i], up[1][i], up[2][i], up[3][i], order, pp);
287 }
288 }
289
290 return pp;
291 }
292
generatePointsPriCGNS(int order,bool complete)293 std::vector<SVector3> generatePointsPriCGNS(int order, bool complete)
294 {
295 std::vector<SVector3> pp;
296
297 if(order == 0) {
298 pp.push_back(SVector3(1. / 3., 1. / 3., 0.));
299 return pp;
300 }
301
302 // principal vertices
303 pp.push_back(SVector3(0, 0, -1));
304 pp.push_back(SVector3(1, 0, -1));
305 pp.push_back(SVector3(0, 1, -1));
306 pp.push_back(SVector3(0, 0, 1));
307 pp.push_back(SVector3(1, 0, 1));
308 pp.push_back(SVector3(0, 1, 1));
309
310 // internal points in edges of base triangle
311 for(int i = 0; i < 3; i++)
312 addEdgePointsCGNS(pp[i], pp[(i + 1) % 3], order, pp);
313
314 // internal points in upstanding edges
315 std::vector<SVector3> edge[3]; // keep for definition of volume pp
316 for(int i = 0; i < 3; i++) {
317 addEdgePointsCGNS(pp[i], pp[i + 3], order, edge[i]);
318 pp.insert(pp.end(), edge[i].begin(), edge[i].end());
319 }
320
321 // internal points in edges of top triangle
322 for(int i = 0; i < 3; i++)
323 addEdgePointsCGNS(pp[i + 3], pp[(i + 1) % 3 + 3], order, pp);
324
325 if(complete) {
326 // internal vertices for base triangle
327 addTriPointsCGNS(pp[0], pp[1], pp[2], order, pp);
328
329 // internal vertices for upstanding quadrilaterals
330 for(int i = 0; i < 3; i++) {
331 addQuaPointsCGNS(pp[i], pp[(i + 1) % 3], pp[(i + 1) % 3 + 3], pp[i + 3],
332 order, pp);
333 }
334
335 // internal points for top triangle
336 addTriPointsCGNS(pp[3], pp[4], pp[5], order, pp);
337
338 // internal points in the volume as a succession of "triangles"
339 for(int o = 0; o < order - 1; o++) {
340 addTriPointsCGNS(edge[0][o], edge[1][o], edge[2][o], order, pp);
341 }
342 }
343
344 return pp;
345 }
346
347 // WARNING: incomplete pyramid order 2 is wrong
generatePointsPyrCGNS(int order,bool complete)348 std::vector<SVector3> generatePointsPyrCGNS(int order, bool complete)
349 {
350 std::vector<SVector3> pp;
351
352 if(order == 0) {
353 pp.push_back(SVector3(0., 0., 0.));
354 return pp;
355 }
356
357 // principal vertices
358 pp.push_back(SVector3(-1, -1, 0));
359 pp.push_back(SVector3(1, -1, 0));
360 pp.push_back(SVector3(1, 1, 0));
361 pp.push_back(SVector3(-1, 1, 0));
362 pp.push_back(SVector3(0, 0, 1));
363
364 // internal points in edges of base quadrilateral
365 for(int i = 0; i < 4; i++)
366 addEdgePointsCGNS(pp[i], pp[(i + 1) % 4], order, pp);
367
368 // internal points in upstanding edges
369 for(int i = 0; i < 4; i++) addEdgePointsCGNS(pp[i], pp[4], order, pp);
370
371 // internal points in base quadrilateral
372 addQuaPointsCGNS(pp[0], pp[1], pp[2], pp[3], order, pp);
373
374 // internal points in upstanding triangles
375 for(int i = 0; i < 4; i++)
376 addTriPointsCGNS(pp[i], pp[(i + 1) % 4], pp[4], order, pp);
377
378 // internal points as an internal pyramid of order p-3
379 if(order > 2) {
380 std::vector<SVector3> pyr = generatePointsPyrCGNS(order - 3, true);
381 SVector3 offset(0, 0, 1. / order);
382 double scale = double(order - 3) / double(order);
383 for(size_t i = 0; i < pyr.size(); ++i)
384 pp.push_back((pyr[i] * scale) + offset);
385 }
386
387 return pp;
388 }
389
generatePointsCGNS(int parentType,int order,bool complete)390 fullMatrix<double> generatePointsCGNS(int parentType, int order,
391 bool complete)
392 {
393 std::vector<SVector3> pts;
394
395 switch(parentType) {
396 case TYPE_LIN: pts = generatePointsEdgeCGNS(order); break;
397 case TYPE_TRI: pts = generatePointsTriCGNS(order, complete); break;
398 case TYPE_QUA: pts = generatePointsQuaCGNS(order, complete); break;
399 case TYPE_TET: pts = generatePointsTetCGNS(order, complete); break;
400 case TYPE_HEX: pts = generatePointsHexCGNS(order, complete); break;
401 case TYPE_PRI: pts = generatePointsPriCGNS(order, complete); break;
402 case TYPE_PYR: pts = generatePointsPyrCGNS(order, complete); break;
403 default:
404 Msg::Error(
405 "%s (%i) : Error CGNS element %s of order %i not yet implemented",
406 __FILE__, __LINE__, ElementType::nameOfParentType(parentType).c_str(),
407 order);
408 }
409
410 size_t dim = 0;
411 switch(parentType) {
412 case TYPE_PNT: dim = 3; break;
413 case TYPE_LIN: dim = 1; break;
414 case TYPE_TRI:
415 case TYPE_QUA: dim = 2; break;
416 case TYPE_TET:
417 case TYPE_HEX:
418 case TYPE_PRI:
419 case TYPE_PYR: dim = 3; break;
420 }
421
422 fullMatrix<double> ptsCGNS(pts.size(), dim);
423 for(size_t i = 0; i < pts.size(); i++) {
424 for(size_t d = 0; d < dim; d++) ptsCGNS(i, d) = pts[i][d];
425 }
426
427 return ptsCGNS;
428 }
429
cgns2MshNodeIndexInit(int mshTag)430 std::vector<int> cgns2MshNodeIndexInit(int mshTag)
431 {
432 const nodalBasis *nb = BasisFactory::getNodalBasis(mshTag);
433 const int nNode = nb->points.size1();
434
435 std::vector<int> nodes(nNode);
436
437 switch(mshTag) {
438 case MSH_PNT:
439 case MSH_LIN_2:
440 case MSH_TRI_3:
441 case MSH_QUA_4:
442 case MSH_TET_4:
443 case MSH_HEX_8:
444 case MSH_PRI_6:
445 case MSH_PYR_5: // linear elements: same ordering between Gmsh and CGNS
446 for(int i = 0; i < nNode; i++) nodes[i] = i;
447 break;
448 default: // high-order elements: get reordering by comparing points
449 int parent = ElementType::getParentType(mshTag);
450 int order = ElementType::getOrder(mshTag);
451 bool complete = (ElementType::getSerendipity(mshTag) <= 1);
452 fullMatrix<double> ptsCGNS = generatePointsCGNS(parent, order, complete);
453 const bool reorderOK = computeReordering(ptsCGNS, nb->points, nodes);
454 if(!reorderOK) {
455 Msg::Error("%s (%i) : Error in computation of reordering between Gmsh "
456 "and CGNS element nodes for MSH type %i",
457 __FILE__, __LINE__, mshTag);
458 }
459 break;
460 }
461
462 return nodes;
463 }
464
465 // ------------- Conversion of element types between CGNS and Gmsh
466 // -------------
467
msh2CgnsEltTypeInit()468 std::vector<CGNS_ENUMT(ElementType_t)> msh2CgnsEltTypeInit()
469 {
470 std::vector<CGNS_ENUMT(ElementType_t)> cgnsType(
471 MSH_MAX_NUM + 1, CGNS_ENUMV(ElementTypeNull));
472
473 // All orders
474 cgnsType[MSH_PNT] = CGNS_ENUMV(NODE);
475
476 // Linear elements
477 cgnsType[MSH_LIN_2] = CGNS_ENUMV(BAR_2);
478 cgnsType[MSH_TRI_3] = CGNS_ENUMV(TRI_3);
479 cgnsType[MSH_QUA_4] = CGNS_ENUMV(QUAD_4);
480 cgnsType[MSH_TET_4] = CGNS_ENUMV(TETRA_4);
481 cgnsType[MSH_PYR_5] = CGNS_ENUMV(PYRA_5);
482 cgnsType[MSH_PRI_6] = CGNS_ENUMV(PENTA_6);
483 cgnsType[MSH_HEX_8] = CGNS_ENUMV(HEXA_8);
484
485 // Quadratic elements
486 cgnsType[MSH_LIN_3] = CGNS_ENUMV(BAR_3);
487 cgnsType[MSH_TRI_6] = CGNS_ENUMV(TRI_6);
488 cgnsType[MSH_QUA_8] = CGNS_ENUMV(QUAD_8);
489 cgnsType[MSH_QUA_9] = CGNS_ENUMV(QUAD_9);
490 cgnsType[MSH_TET_10] = CGNS_ENUMV(TETRA_10);
491 cgnsType[MSH_PYR_13] = CGNS_ENUMV(PYRA_13);
492 cgnsType[MSH_PYR_14] = CGNS_ENUMV(PYRA_14);
493 cgnsType[MSH_PRI_15] = CGNS_ENUMV(PENTA_15);
494 cgnsType[MSH_PRI_18] = CGNS_ENUMV(PENTA_18);
495 cgnsType[MSH_HEX_20] = CGNS_ENUMV(HEXA_20);
496 cgnsType[MSH_HEX_27] = CGNS_ENUMV(HEXA_27);
497
498 // Cubic elements
499 cgnsType[MSH_LIN_4] = CGNS_ENUMV(BAR_4);
500 cgnsType[MSH_TRI_9] = CGNS_ENUMV(TRI_9);
501 cgnsType[MSH_TRI_10] = CGNS_ENUMV(TRI_10);
502 cgnsType[MSH_QUA_12] = CGNS_ENUMV(QUAD_12);
503 cgnsType[MSH_QUA_16] = CGNS_ENUMV(QUAD_16);
504 cgnsType[MSH_TET_16] = CGNS_ENUMV(TETRA_16);
505 cgnsType[MSH_TET_20] = CGNS_ENUMV(TETRA_20);
506 cgnsType[MSH_PYR_21] = CGNS_ENUMV(PYRA_21);
507 cgnsType[MSH_PYR_29] = CGNS_ENUMV(PYRA_29);
508 cgnsType[MSH_PYR_30] = CGNS_ENUMV(PYRA_30);
509 cgnsType[MSH_PRI_24] = CGNS_ENUMV(PENTA_24);
510 // cgnsType[MSH_PRI_38] = CGNS_ENUMV(PENTA_38);
511 cgnsType[MSH_PRI_40] = CGNS_ENUMV(PENTA_40);
512 cgnsType[MSH_HEX_32] = CGNS_ENUMV(HEXA_32);
513 cgnsType[MSH_HEX_56] = CGNS_ENUMV(HEXA_56);
514 cgnsType[MSH_HEX_64] = CGNS_ENUMV(HEXA_64);
515
516 // Quartic elements
517 cgnsType[MSH_LIN_5] = CGNS_ENUMV(BAR_5);
518 cgnsType[MSH_TRI_12] = CGNS_ENUMV(TRI_12);
519 cgnsType[MSH_TRI_15] = CGNS_ENUMV(TRI_15);
520 cgnsType[MSH_QUA_16] = CGNS_ENUMV(QUAD_16);
521 cgnsType[MSH_QUA_25] = CGNS_ENUMV(QUAD_25);
522 cgnsType[MSH_TET_22] = CGNS_ENUMV(TETRA_22);
523 cgnsType[MSH_TET_34] = CGNS_ENUMV(TETRA_34);
524 cgnsType[MSH_TET_35] = CGNS_ENUMV(TETRA_35);
525 cgnsType[MSH_PYR_29] = CGNS_ENUMV(PYRA_29);
526 // cgnsType[MSH_PYR_50] = CGNS_ENUMV(PYRA_50);
527 cgnsType[MSH_PYR_55] = CGNS_ENUMV(PYRA_55);
528 cgnsType[MSH_PRI_33] = CGNS_ENUMV(PENTA_33);
529 // cgnsType[MSH_PRI_66] = CGNS_ENUMV(PENTA_66);
530 cgnsType[MSH_PRI_75] = CGNS_ENUMV(PENTA_75);
531 cgnsType[MSH_HEX_44] = CGNS_ENUMV(HEXA_44);
532 // cgnsType[MSH_HEX_98] = CGNS_ENUMV(HEXA_98);
533 cgnsType[MSH_HEX_125] = CGNS_ENUMV(HEXA_125);
534
535 return cgnsType;
536 }
537
cgns2MshEltTypeInit()538 std::vector<int> cgns2MshEltTypeInit()
539 {
540 std::vector<int> mshType(NofValidElementTypes, 0);
541
542 // All orders
543 mshType[CGNS_ENUMV(NODE)] = MSH_PNT;
544
545 // Linear elements
546 mshType[CGNS_ENUMV(BAR_2)] = MSH_LIN_2;
547 mshType[CGNS_ENUMV(TRI_3)] = MSH_TRI_3;
548 mshType[CGNS_ENUMV(QUAD_4)] = MSH_QUA_4;
549 mshType[CGNS_ENUMV(TETRA_4)] = MSH_TET_4;
550 mshType[CGNS_ENUMV(PYRA_5)] = MSH_PYR_5;
551 mshType[CGNS_ENUMV(PENTA_6)] = MSH_PRI_6;
552 mshType[CGNS_ENUMV(HEXA_8)] = MSH_HEX_8;
553
554 // Quadratic elements
555 mshType[CGNS_ENUMV(BAR_3)] = MSH_LIN_3;
556 mshType[CGNS_ENUMV(TRI_6)] = MSH_TRI_6;
557 mshType[CGNS_ENUMV(QUAD_8)] = MSH_QUA_8;
558 mshType[CGNS_ENUMV(QUAD_9)] = MSH_QUA_9;
559 mshType[CGNS_ENUMV(TETRA_10)] = MSH_TET_10;
560 mshType[CGNS_ENUMV(PYRA_13)] = MSH_PYR_13;
561 mshType[CGNS_ENUMV(PYRA_14)] = MSH_PYR_14;
562 mshType[CGNS_ENUMV(PENTA_15)] = MSH_PRI_15;
563 mshType[CGNS_ENUMV(PENTA_18)] = MSH_PRI_18;
564 mshType[CGNS_ENUMV(HEXA_20)] = MSH_HEX_20;
565 mshType[CGNS_ENUMV(HEXA_27)] = MSH_HEX_27;
566
567 // Cubic elements
568 mshType[CGNS_ENUMV(BAR_4)] = MSH_LIN_4;
569 mshType[CGNS_ENUMV(TRI_9)] = MSH_TRI_9;
570 mshType[CGNS_ENUMV(TRI_10)] = MSH_TRI_10;
571 mshType[CGNS_ENUMV(QUAD_12)] = MSH_QUA_12;
572 mshType[CGNS_ENUMV(QUAD_16)] = MSH_QUA_16;
573 mshType[CGNS_ENUMV(TETRA_16)] = MSH_TET_16;
574 mshType[CGNS_ENUMV(TETRA_20)] = MSH_TET_20;
575 mshType[CGNS_ENUMV(PYRA_21)] = MSH_PYR_21;
576 mshType[CGNS_ENUMV(PYRA_29)] = MSH_PYR_29;
577 mshType[CGNS_ENUMV(PYRA_30)] = MSH_PYR_30;
578 mshType[CGNS_ENUMV(PENTA_24)] = MSH_PRI_24;
579 // mshType[CGNS_ENUMV(PENTA_38)] = MSH_PRI_38;
580 mshType[CGNS_ENUMV(PENTA_40)] = MSH_PRI_40;
581 mshType[CGNS_ENUMV(HEXA_32)] = MSH_HEX_32;
582 mshType[CGNS_ENUMV(HEXA_56)] = MSH_HEX_56;
583 mshType[CGNS_ENUMV(HEXA_64)] = MSH_HEX_64;
584
585 // Quartic elements
586 mshType[CGNS_ENUMV(BAR_5)] = MSH_LIN_5;
587 mshType[CGNS_ENUMV(TRI_12)] = MSH_TRI_12;
588 mshType[CGNS_ENUMV(TRI_15)] = MSH_TRI_15;
589 mshType[CGNS_ENUMV(QUAD_16)] = MSH_QUA_16;
590 mshType[CGNS_ENUMV(QUAD_25)] = MSH_QUA_25;
591 mshType[CGNS_ENUMV(TETRA_22)] = MSH_TET_22;
592 mshType[CGNS_ENUMV(TETRA_34)] = MSH_TET_34;
593 mshType[CGNS_ENUMV(TETRA_35)] = MSH_TET_35;
594 mshType[CGNS_ENUMV(PYRA_29)] = MSH_PYR_29;
595 // mshType[CGNS_ENUMV(PYRA_50)] = MSH_PYR_50;
596 mshType[CGNS_ENUMV(PYRA_55)] = MSH_PYR_55;
597 mshType[CGNS_ENUMV(PENTA_33)] = MSH_PRI_33;
598 // mshType[CGNS_ENUMV(PENTA_66)] = MSH_PRI_66;
599 mshType[CGNS_ENUMV(PENTA_75)] = MSH_PRI_75;
600 mshType[CGNS_ENUMV(HEXA_44)] = MSH_HEX_44;
601 // mshType[CGNS_ENUMV(HEXA_98)] = MSH_HEX_98;
602 mshType[CGNS_ENUMV(HEXA_125)] = MSH_HEX_125;
603
604 return mshType;
605 }
606
607 } // namespace
608
609 // msh to CGNS element type
msh2CgnsEltType(int mshTag)610 CGNS_ENUMT(ElementType_t) msh2CgnsEltType(int mshTag)
611 {
612 static std::vector<CGNS_ENUMT(ElementType_t)> cgnsType =
613 msh2CgnsEltTypeInit();
614
615 if(mshTag >= static_cast<int>(cgnsType.size()))
616 return CGNS_ENUMV(ElementTypeNull);
617 return cgnsType[mshTag];
618 }
619
620 // CGNS to msh element type
cgns2MshEltType(CGNS_ENUMT (ElementType_t)cgnsType)621 int cgns2MshEltType(CGNS_ENUMT(ElementType_t) cgnsType)
622 {
623 static std::vector<int> mshType = cgns2MshEltTypeInit();
624
625 if(cgnsType >= static_cast<int>(mshType.size())) return 0;
626 return mshType[cgnsType];
627 }
628
629 // CGNS to msh node ordering
cgns2MshNodeIndex(int mshTag)630 std::vector<int> &cgns2MshNodeIndex(int mshTag)
631 {
632 static std::vector<std::vector<int> > mshInd(MSH_MAX_NUM + 1);
633 static std::vector<int> dumInd;
634
635 if(mshTag > MSH_MAX_NUM) return dumInd;
636
637 if(mshInd[mshTag].size() == 0) {
638 mshInd[mshTag] = cgns2MshNodeIndexInit(mshTag);
639 }
640
641 return mshInd[mshTag];
642 }
643
644 // // msh to CGNS node ordering
645 // std::vector<int> &msh2CgnsNodeIndex(int mshTag)
646 // {
647 // static std::vector<std::vector<int> > cgnsInd(MSH_MAX_NUM+1);
648 // static std::vector<int> dumInd;
649
650 // if(mshTag > MSH_MAX_NUM) return dumInd;
651
652 // if(cgnsInd[mshTag].size() == 0) {
653 // std::vector<int> &mshInd = cgns2MshNodeIndex(mshTag);
654 // std::vector<int> &ind = cgnsInd[mshTag];
655 // for(int iV = 0; iV < mshInd.size(); iV++) ind[mshInd[iV]] = iV;
656 // }
657
658 // return cgnsInd[mshTag];
659 // }
660
msh2CgnsReferenceElement(int mshType,const fullMatrix<double> & mshPts,std::vector<double> & u,std::vector<double> & v,std::vector<double> & w)661 void msh2CgnsReferenceElement(int mshType, const fullMatrix<double> &mshPts,
662 std::vector<double> &u, std::vector<double> &v,
663 std::vector<double> &w)
664 {
665 int parentType = ElementType::getParentType(mshType);
666
667 switch(parentType) {
668 case TYPE_PNT: u[0] = mshPts(0, 0); break;
669 case TYPE_LIN: // Gmsh and CGNS both in [-1, 1]
670 for(int i = 0; i < mshPts.size1(); i++) { u[i] = mshPts(i, 0); }
671 break;
672 case TYPE_QUA: // Gmsh and CGNS both in [-1, 1]
673 for(int i = 0; i < mshPts.size1(); i++) {
674 u[i] = mshPts(i, 0);
675 v[i] = mshPts(i, 1);
676 }
677 break;
678 case TYPE_HEX: // Gmsh and CGNS both in [-1, 1]
679 for(int i = 0; i < mshPts.size1(); i++) {
680 u[i] = mshPts(i, 0);
681 v[i] = mshPts(i, 1);
682 w[i] = mshPts(i, 2);
683 }
684 break;
685 case TYPE_TRI: // Gmsh in [0, 1], CGNS in [-1, 1]
686 for(int i = 0; i < mshPts.size1(); i++) {
687 u[i] = 2. * mshPts(i, 0) - 1.;
688 v[i] = 2. * mshPts(i, 1) - 1.;
689 }
690 break;
691 case TYPE_TET: // Gmsh in [0, 1], CGNS in [-1, 1]
692 for(int i = 0; i < mshPts.size1(); i++) {
693 u[i] = 2. * mshPts(i, 0) - 1.;
694 v[i] = 2. * mshPts(i, 1) - 1.;
695 w[i] = 2. * mshPts(i, 2) - 1.;
696 }
697 break;
698 case TYPE_PRI: // uv: Gmsh in [0, 1] and CGNS in [-1, 1], w: both in [-1, 1]
699 for(int i = 0; i < mshPts.size1(); i++) {
700 u[i] = 2. * mshPts(i, 0) - 1.;
701 v[i] = 2. * mshPts(i, 1) - 1.;
702 w[i] = mshPts(i, 2);
703 }
704 break;
705 case TYPE_PYR: // uv: both in [-1, 1], w: Gmsh in [0, 1] and CGNS in [-1, 1]
706 for(int i = 0; i < mshPts.size1(); i++) {
707 u[i] = mshPts(i, 0);
708 v[i] = mshPts(i, 1);
709 w[i] = 2. * mshPts(i, 2) - 1.;
710 }
711 break;
712 default:
713 Msg::Error("%s (%i) : Error CGNS element %s not yet implemented", __FILE__,
714 __LINE__, ElementType::nameOfParentType(parentType).c_str());
715 }
716 }
717
computeReordering(fullMatrix<double> src,fullMatrix<double> dest,std::vector<int> & ind)718 bool computeReordering(fullMatrix<double> src, fullMatrix<double> dest,
719 std::vector<int> &ind)
720 {
721 static const double TOLSQ = 1e-10 * 1e-10;
722
723 const size_t nNode = src.size1(), dim = src.size2();
724 ind.resize(nNode, -1);
725
726 // Loop over src and dest nodes and check distance
727 std::set<int> found;
728 for(size_t iSrc = 0; iSrc < nNode; iSrc++) {
729 for(size_t iDest = 0; iDest < nNode; iDest++) {
730 const double xx = src(iSrc, 0) - dest(iDest, 0);
731 const double yy = (dim > 1) ? src(iSrc, 1) - dest(iDest, 1) : 0.;
732 const double zz = (dim > 2) ? src(iSrc, 2) - dest(iDest, 2) : 0.;
733 const double dd = xx * xx + yy * yy + zz * zz;
734 if(dd < TOLSQ) {
735 ind[iSrc] = iDest;
736 found.insert(iDest);
737 break;
738 }
739 }
740 }
741
742 // Check bijectivity
743 return (found.size() == nNode);
744 }
745
cgnsString(const std::string & s,std::string::size_type maxLength)746 std::string cgnsString(const std::string &s, std::string::size_type maxLength)
747 {
748 std::string s2(s);
749 if(s2.size() > maxLength) s2.resize(maxLength);
750 return s2;
751 }
752
753 #endif // HAVE_LIBCGNS
754