Actual source code: plexgmsh.c
petsc-3.10.5 2019-03-28
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
4: /*@C
5: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
7: + comm - The MPI communicator
8: . filename - Name of the Gmsh file
9: - interpolate - Create faces and edges in the mesh
11: Output Parameter:
12: . dm - The DM object representing the mesh
14: Level: beginner
16: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
17: @*/
18: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
19: {
20: PetscViewer viewer;
21: PetscMPIInt rank;
22: int fileType;
23: PetscViewerType vtype;
24: PetscErrorCode ierr;
27: MPI_Comm_rank(comm, &rank);
29: /* Determine Gmsh file type (ASCII or binary) from file header */
30: if (!rank) {
31: PetscViewer vheader;
32: char line[PETSC_MAX_PATH_LEN];
33: PetscBool match;
34: int snum;
35: float version;
37: PetscViewerCreate(PETSC_COMM_SELF, &vheader);
38: PetscViewerSetType(vheader, PETSCVIEWERASCII);
39: PetscViewerFileSetMode(vheader, FILE_MODE_READ);
40: PetscViewerFileSetName(vheader, filename);
41: /* Read only the first two lines of the Gmsh file */
42: PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);
43: PetscStrncmp(line, "$MeshFormat", sizeof(line), &match);
44: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
45: PetscViewerRead(vheader, line, 2, NULL, PETSC_STRING);
46: snum = sscanf(line, "%f %d", &version, &fileType);
47: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
48: if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
49: PetscViewerDestroy(&vheader);
50: }
51: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
52: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
54: /* Create appropriate viewer and build plex */
55: PetscViewerCreate(comm, &viewer);
56: PetscViewerSetType(viewer, vtype);
57: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
58: PetscViewerFileSetName(viewer, filename);
59: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
60: PetscViewerDestroy(&viewer);
61: return(0);
62: }
64: static PetscErrorCode DMPlexCreateGmsh_ReadNodes(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, int shift, int numVertices, double **coordinates)
65: {
66: int v,nid;
70: PetscMalloc1(numVertices*3, coordinates);
71: for (v = 0; v < numVertices; ++v) {
72: double *xyz = *coordinates + v*3;
73: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
74: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
75: if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift);
76: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
77: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
78: }
79: return(0);
80: }
82: static PetscErrorCode DMPlexCreateGmsh_ReadElements(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, PETSC_UNUSED int shift, int numCells, GmshElement **gmsh_elems)
83: {
84: GmshElement *elements;
85: int i, c, p, ibuf[1+4+512];
86: int cellType, dim, numNodes, numNodesIgnore, numElem, numTags;
90: PetscMalloc1(numCells, &elements);
91: for (c = 0; c < numCells;) {
92: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
93: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
94: if (binary) {
95: cellType = ibuf[0];
96: numElem = ibuf[1];
97: numTags = ibuf[2];
98: } else {
99: elements[c].id = ibuf[0];
100: cellType = ibuf[1];
101: numTags = ibuf[2];
102: numElem = 1;
103: }
104: /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */
105: numNodesIgnore = 0;
106: switch (cellType) {
107: case 1: /* 2-node line */
108: dim = 1;
109: numNodes = 2;
110: break;
111: case 2: /* 3-node triangle */
112: dim = 2;
113: numNodes = 3;
114: break;
115: case 3: /* 4-node quadrangle */
116: dim = 2;
117: numNodes = 4;
118: break;
119: case 4: /* 4-node tetrahedron */
120: dim = 3;
121: numNodes = 4;
122: break;
123: case 5: /* 8-node hexahedron */
124: dim = 3;
125: numNodes = 8;
126: break;
127: case 6: /* 6-node wedge */
128: dim = 3;
129: numNodes = 6;
130: break;
131: case 8: /* 3-node 2nd order line */
132: dim = 1;
133: numNodes = 2;
134: numNodesIgnore = 1;
135: break;
136: case 9: /* 6-node 2nd order triangle */
137: dim = 2;
138: numNodes = 3;
139: numNodesIgnore = 3;
140: break;
141: case 13: /* 18-node 2nd wedge */
142: dim = 3;
143: numNodes = 6;
144: numNodesIgnore = 12;
145: break;
146: case 15: /* 1-node vertex */
147: dim = 0;
148: numNodes = 1;
149: break;
150: case 7: /* 5-node pyramid */
151: case 10: /* 9-node 2nd order quadrangle */
152: case 11: /* 10-node 2nd order tetrahedron */
153: case 12: /* 27-node 2nd order hexhedron */
154: case 14: /* 14-node 2nd order pyramid */
155: default:
156: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
157: }
158: if (binary) {
159: const int nint = 1 + numTags + numNodes + numNodesIgnore;
160: /* Loop over element blocks */
161: for (i = 0; i < numElem; ++i, ++c) {
162: PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
163: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, nint);}
164: elements[c].dim = dim;
165: elements[c].numNodes = numNodes;
166: elements[c].numTags = numTags;
167: elements[c].id = ibuf[0];
168: elements[c].cellType = cellType;
169: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
170: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
171: }
172: } else {
173: elements[c].dim = dim;
174: elements[c].numNodes = numNodes;
175: elements[c].numTags = numTags;
176: elements[c].cellType = cellType;
177: PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
178: PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
179: PetscViewerRead(viewer, ibuf, numNodesIgnore, NULL, PETSC_ENUM);
180: c++;
181: }
182: }
183: *gmsh_elems = elements;
184: return(0);
185: }
187: /*@
188: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
190: Collective on comm
192: Input Parameters:
193: + comm - The MPI communicator
194: . viewer - The Viewer associated with a Gmsh file
195: - interpolate - Create faces and edges in the mesh
197: Output Parameter:
198: . dm - The DM object representing the mesh
200: Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format
201: and http://www.geuz.org/gmsh/doc/texinfo/#MSH-binary-file-format
203: Level: beginner
205: .keywords: mesh,Gmsh
206: .seealso: DMPLEX, DMCreate()
207: @*/
208: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
209: {
210: PetscViewer parentviewer = NULL;
211: double *coordsIn = NULL;
212: GmshElement *gmsh_elem = NULL;
213: PetscSection coordSection;
214: Vec coordinates;
215: PetscBT periodicV = NULL, periodicC = NULL;
216: PetscScalar *coords;
217: PetscInt dim = 0, embedDim = -1, coordSize, c, v, d, r, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL, cMax = PETSC_DETERMINE;
218: int i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum, shift = 1;
219: PetscMPIInt rank;
220: char line[PETSC_MAX_PATH_LEN];
221: PetscBool binary, byteSwap = PETSC_FALSE, zerobase = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE;
222: PetscBool enable_hybrid = PETSC_FALSE;
226: MPI_Comm_rank(comm, &rank);
227: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_hybrid", &enable_hybrid, NULL);
228: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_periodic", &periodic, NULL);
229: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_use_marker", &usemarker, NULL);
230: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_zero_base", &zerobase, NULL);
231: PetscOptionsGetInt(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_spacedim", &embedDim, NULL);
232: if (zerobase) shift = 0;
234: DMCreate(comm, dm);
235: DMSetType(*dm, DMPLEX);
236: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
238: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
240: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
241: if (binary) {
242: parentviewer = viewer;
243: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
244: }
246: if (!rank) {
247: PetscBool match, hybrid;
248: int fileType, dataSize;
249: float version;
251: /* Read header */
252: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
253: PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
254: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
255: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
256: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
257: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
258: if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
259: if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
260: if (binary) {
261: int checkInt;
262: PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
263: if (checkInt != 1) {
264: PetscByteSwap(&checkInt, PETSC_ENUM, 1);
265: if (checkInt == 1) byteSwap = PETSC_TRUE;
266: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
267: }
268: } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
269: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
270: PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
271: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
273: /* OPTIONAL Read physical names */
274: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
275: PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);
276: if (match) {
277: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
278: snum = sscanf(line, "%d", &numRegions);
279: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
280: for (r = 0; r < numRegions; ++r) {
281: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
282: }
283: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
284: PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);
285: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
286: /* Initial read for vertex section */
287: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
288: }
290: /* Read vertices */
291: PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
292: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
293: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
294: snum = sscanf(line, "%d", &numVertices);
295: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
296: DMPlexCreateGmsh_ReadNodes(viewer, binary, byteSwap, shift, numVertices, &coordsIn);
297: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
298: PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
299: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
301: /* Read cells */
302: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
303: PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
304: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
305: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
306: snum = sscanf(line, "%d", &numCells);
307: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
308: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
309: file contents multiple times to figure out the true number of cells and facets
310: in the given mesh. To make this more efficient we read the file contents only
311: once and store them in memory, while determining the true number of cells. */
312: DMPlexCreateGmsh_ReadElements(viewer, binary, byteSwap, shift, numCells, &gmsh_elem);
313: hybrid = PETSC_FALSE;
314: for (trueNumCells = 0, c = 0; c < numCells; ++c) {
315: int on = -1;
316: if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
317: if (gmsh_elem[c].dim == dim) {hybrid = (trueNumCells ? (on != gmsh_elem[c].cellType ? on = gmsh_elem[c].cellType,PETSC_TRUE : hybrid) : (on = gmsh_elem[c].cellType, PETSC_FALSE) ); trueNumCells++;}
318: /* wedges always indicate an hybrid mesh in PLEX */
319: if (on == 6 || on == 13) hybrid = PETSC_TRUE;
320: }
321: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
322: PetscStrncmp(line, "$EndElements", 12, &match);
323: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
325: /* Renumber cells for hybrid grids */
326: if (hybrid && enable_hybrid) {
327: PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
328: PetscInt cell, tn, *tp;
329: int n1 = 0,n2 = 0;
331: PetscMalloc1(trueNumCells, &hybridCells1);
332: PetscMalloc1(trueNumCells, &hybridCells2);
333: for (cell = 0, c = 0; c < numCells; ++c) {
334: if (gmsh_elem[c].dim == dim) {
335: if (!n1) n1 = gmsh_elem[c].cellType;
336: else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;
338: if (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
339: else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
340: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
341: cell++;
342: }
343: }
345: switch (n1) {
346: case 2: /* triangles */
347: case 9:
348: switch (n2) {
349: case 0: /* single type mesh */
350: case 3: /* quads */
351: case 10:
352: break;
353: default:
354: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
355: }
356: break;
357: case 3:
358: case 10:
359: switch (n2) {
360: case 0: /* single type mesh */
361: case 2: /* swap since we list simplices first */
362: case 9:
363: tn = hc1;
364: hc1 = hc2;
365: hc2 = tn;
367: tp = hybridCells1;
368: hybridCells1 = hybridCells2;
369: hybridCells2 = tp;
370: break;
371: default:
372: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
373: }
374: break;
375: case 4: /* tetrahedra */
376: case 11:
377: switch (n2) {
378: case 0: /* single type mesh */
379: case 6: /* wedges */
380: case 13:
381: break;
382: default:
383: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
384: }
385: break;
386: case 6:
387: case 13:
388: switch (n2) {
389: case 0: /* single type mesh */
390: case 4: /* swap since we list simplices first */
391: case 11:
392: tn = hc1;
393: hc1 = hc2;
394: hc2 = tn;
396: tp = hybridCells1;
397: hybridCells1 = hybridCells2;
398: hybridCells2 = tp;
399: break;
400: default:
401: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
402: }
403: break;
404: default:
405: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
406: }
407: cMax = hc1;
408: PetscMalloc1(trueNumCells, &hybridMap);
409: for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
410: for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
411: PetscFree(hybridCells1);
412: PetscFree(hybridCells2);
413: }
415: /* OPTIONAL Read periodic section */
416: if (periodic) {
417: PetscInt pVert, *periodicMapT, *aux;
418: int numPeriodic;
420: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
421: PetscStrncmp(line, "$Periodic", 9, &match);
422: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
423: PetscMalloc1(numVertices, &periodicMapT);
424: PetscBTCreate(numVertices, &periodicV);
425: for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
426: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
427: snum = sscanf(line, "%d", &numPeriodic);
428: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
429: for (i = 0; i < numPeriodic; i++) {
430: int j, edim = -1, slaveTag = -1, masterTag = -1, nNodes, slaveNode, masterNode;
432: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
433: snum = sscanf(line, "%d %d %d", &edim, &slaveTag, &masterTag); /* slaveTag and masterTag are unused */
434: if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
435: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
436: snum = sscanf(line, "%d", &nNodes);
437: if (snum != 1) { /* discard tranformation and try again */
438: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
439: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
440: snum = sscanf(line, "%d", &nNodes);
441: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
442: }
443: for (j = 0; j < nNodes; j++) {
444: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
445: snum = sscanf(line, "%d %d", &slaveNode, &masterNode);
446: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
447: periodicMapT[slaveNode - shift] = masterNode - shift;
448: PetscBTSet(periodicV, slaveNode - shift);
449: PetscBTSet(periodicV, masterNode - shift);
450: }
451: }
452: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
453: PetscStrncmp(line, "$EndPeriodic", 12, &match);
454: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
455: /* we may have slaves of slaves */
456: for (i = 0; i < numVertices; i++) {
457: while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
458: periodicMapT[i] = periodicMapT[periodicMapT[i]];
459: }
460: }
461: /* periodicMap : from old to new numbering (periodic vertices excluded)
462: periodicMapI: from new to old numbering */
463: PetscMalloc1(numVertices, &periodicMap);
464: PetscMalloc1(numVertices, &periodicMapI);
465: PetscMalloc1(numVertices, &aux);
466: for (i = 0, pVert = 0; i < numVertices; i++) {
467: if (periodicMapT[i] != i) {
468: pVert++;
469: } else {
470: aux[i] = i - pVert;
471: periodicMapI[i - pVert] = i;
472: }
473: }
474: for (i = 0 ; i < numVertices; i++) {
475: periodicMap[i] = aux[periodicMapT[i]];
476: }
477: PetscFree(periodicMapT);
478: PetscFree(aux);
479: /* remove periodic vertices */
480: numVertices = numVertices - pVert;
481: }
482: }
484: if (parentviewer) {
485: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
486: }
488: /* Allocate the cell-vertex mesh */
489: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
490: for (cell = 0, c = 0; c < numCells; ++c) {
491: if (gmsh_elem[c].dim == dim) {
492: DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);
493: cell++;
494: }
495: }
496: DMSetUp(*dm);
497: /* Add cell-vertex connections */
498: for (cell = 0, c = 0; c < numCells; ++c) {
499: if (gmsh_elem[c].dim == dim) {
500: PetscInt pcone[8], corner;
501: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
502: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
503: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
504: }
505: if (dim == 3) {
506: /* Tetrahedra are inverted */
507: if (gmsh_elem[c].cellType == 4) {
508: PetscInt tmp = pcone[0];
509: pcone[0] = pcone[1];
510: pcone[1] = tmp;
511: }
512: /* Hexahedra are inverted */
513: if (gmsh_elem[c].cellType == 5) {
514: PetscInt tmp = pcone[1];
515: pcone[1] = pcone[3];
516: pcone[3] = tmp;
517: }
518: /* Prisms are inverted */
519: if (gmsh_elem[c].cellType == 6) {
520: PetscInt tmp;
522: tmp = pcone[1];
523: pcone[1] = pcone[2];
524: pcone[2] = tmp;
525: tmp = pcone[4];
526: pcone[4] = pcone[5];
527: pcone[5] = tmp;
528: }
529: } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
530: /* quads are input to PLEX as prisms */
531: if (gmsh_elem[c].cellType == 3) {
532: PetscInt tmp = pcone[2];
533: pcone[2] = pcone[3];
534: pcone[3] = tmp;
535: }
536: }
537: DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);
538: cell++;
539: }
540: }
541: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
542: DMSetDimension(*dm, dim);
543: DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
544: DMPlexSymmetrize(*dm);
545: DMPlexStratify(*dm);
546: if (interpolate) {
547: DM idm;
549: DMPlexInterpolate(*dm, &idm);
550: DMDestroy(dm);
551: *dm = idm;
552: }
554: if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
555: if (!rank && usemarker) {
556: PetscInt f, fStart, fEnd;
558: DMCreateLabel(*dm, "marker");
559: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
560: for (f = fStart; f < fEnd; ++f) {
561: PetscInt suppSize;
563: DMPlexGetSupportSize(*dm, f, &suppSize);
564: if (suppSize == 1) {
565: PetscInt *cone = NULL, coneSize, p;
567: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
568: for (p = 0; p < coneSize; p += 2) {
569: DMSetLabelValue(*dm, "marker", cone[p], 1);
570: }
571: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
572: }
573: }
574: }
576: if (!rank) {
577: PetscInt vStart, vEnd;
579: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
580: for (cell = 0, c = 0; c < numCells; ++c) {
582: /* Create face sets */
583: if (interpolate && gmsh_elem[c].dim == dim-1) {
584: const PetscInt *join;
585: PetscInt joinSize, pcone[8], corner;
586: /* Find the relevant facet with vertex joins */
587: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
588: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
589: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
590: }
591: DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
592: if (joinSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for GMsh element %d (Plex cell %D)", gmsh_elem[c].id, c);
593: DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
594: DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
595: }
597: /* Create cell sets */
598: if (gmsh_elem[c].dim == dim) {
599: if (gmsh_elem[c].numTags > 0) {
600: DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
601: }
602: cell++;
603: }
605: /* Create vertex sets */
606: if (gmsh_elem[c].dim == 0) {
607: if (gmsh_elem[c].numTags > 0) {
608: const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
609: const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
610: DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
611: }
612: }
613: }
614: }
616: /* Create coordinates */
617: if (embedDim < 0) embedDim = dim;
618: DMSetCoordinateDim(*dm, embedDim);
619: DMGetCoordinateSection(*dm, &coordSection);
620: PetscSectionSetNumFields(coordSection, 1);
621: PetscSectionSetFieldComponents(coordSection, 0, embedDim);
622: if (periodicMap) { /* we need to localize coordinates on cells */
623: PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
624: } else {
625: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
626: }
627: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
628: PetscSectionSetDof(coordSection, v, embedDim);
629: PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
630: }
631: if (periodicMap) {
632: PetscBTCreate(trueNumCells, &periodicC);
633: for (cell = 0, c = 0; c < numCells; ++c) {
634: if (gmsh_elem[c].dim == dim) {
635: PetscInt corner;
636: PetscBool pc = PETSC_FALSE;
637: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
638: pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
639: }
640: if (pc) {
641: PetscInt dof = gmsh_elem[c].numNodes*embedDim;
642: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
643: PetscBTSet(periodicC, ucell);
644: PetscSectionSetDof(coordSection, ucell, dof);
645: PetscSectionSetFieldDof(coordSection, ucell, 0, dof);
646: }
647: cell++;
648: }
649: }
650: }
651: PetscSectionSetUp(coordSection);
652: PetscSectionGetStorageSize(coordSection, &coordSize);
653: VecCreate(PETSC_COMM_SELF, &coordinates);
654: PetscObjectSetName((PetscObject) coordinates, "coordinates");
655: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
656: VecSetBlockSize(coordinates, embedDim);
657: VecSetType(coordinates, VECSTANDARD);
658: VecGetArray(coordinates, &coords);
659: if (periodicMap) {
660: PetscInt off;
662: for (cell = 0, c = 0; c < numCells; ++c) {
663: PetscInt pcone[8], corner;
664: if (gmsh_elem[c].dim == dim) {
665: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
666: if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
667: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
668: pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
669: }
670: if (dim == 3) {
671: /* Tetrahedra are inverted */
672: if (gmsh_elem[c].cellType == 4) {
673: PetscInt tmp = pcone[0];
674: pcone[0] = pcone[1];
675: pcone[1] = tmp;
676: }
677: /* Hexahedra are inverted */
678: if (gmsh_elem[c].cellType == 5) {
679: PetscInt tmp = pcone[1];
680: pcone[1] = pcone[3];
681: pcone[3] = tmp;
682: }
683: /* Prisms are inverted */
684: if (gmsh_elem[c].cellType == 6) {
685: PetscInt tmp;
687: tmp = pcone[1];
688: pcone[1] = pcone[2];
689: pcone[2] = tmp;
690: tmp = pcone[4];
691: pcone[4] = pcone[5];
692: pcone[5] = tmp;
693: }
694: } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
695: /* quads are input to PLEX as prisms */
696: if (gmsh_elem[c].cellType == 3) {
697: PetscInt tmp = pcone[2];
698: pcone[2] = pcone[3];
699: pcone[3] = tmp;
700: }
701: }
702: PetscSectionGetOffset(coordSection, ucell, &off);
703: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
704: v = pcone[corner];
705: for (d = 0; d < embedDim; ++d) {
706: coords[off++] = (PetscReal) coordsIn[v*3+d];
707: }
708: }
709: }
710: cell++;
711: }
712: }
713: for (v = 0; v < numVertices; ++v) {
714: PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
715: for (d = 0; d < embedDim; ++d) {
716: coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
717: }
718: }
719: } else {
720: for (v = 0; v < numVertices; ++v) {
721: for (d = 0; d < embedDim; ++d) {
722: coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
723: }
724: }
725: }
726: VecRestoreArray(coordinates, &coords);
727: DMSetCoordinatesLocal(*dm, coordinates);
728: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
730: PetscFree(coordsIn);
731: PetscFree(gmsh_elem);
732: PetscFree(hybridMap);
733: PetscFree(periodicMap);
734: PetscFree(periodicMapI);
735: PetscBTDestroy(&periodicV);
736: PetscBTDestroy(&periodicC);
737: VecDestroy(&coordinates);
739: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
740: return(0);
741: }