Actual source code: plexgmsh.c
petsc-3.8.4 2018-03-24
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, vheader;
21: PetscMPIInt rank;
22: PetscViewerType vtype;
23: char line[PETSC_MAX_PATH_LEN];
24: int snum;
25: PetscBool match;
26: int fT;
27: PetscInt fileType;
28: float version;
29: PetscErrorCode ierr;
32: MPI_Comm_rank(comm, &rank);
33: /* Determine Gmsh file type (ASCII or binary) from file header */
34: PetscViewerCreate(comm, &vheader);
35: PetscViewerSetType(vheader, PETSCVIEWERASCII);
36: PetscViewerFileSetMode(vheader, FILE_MODE_READ);
37: PetscViewerFileSetName(vheader, filename);
38: if (!rank) {
39: /* Read only the first two lines of the Gmsh file */
40: PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);
41: PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
42: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
43: PetscViewerRead(vheader, line, 2, NULL, PETSC_STRING);
44: snum = sscanf(line, "%f %d", &version, &fT);
45: fileType = (PetscInt) fT;
46: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
47: if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
48: }
49: MPI_Bcast(&fileType, 1, MPIU_INT, 0, comm);
50: /* Create appropriate viewer and build plex */
51: if (fileType == 0) vtype = PETSCVIEWERASCII;
52: else vtype = PETSCVIEWERBINARY;
53: PetscViewerCreate(comm, &viewer);
54: PetscViewerSetType(viewer, vtype);
55: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
56: PetscViewerFileSetName(viewer, filename);
57: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
58: PetscViewerDestroy(&viewer);
59: PetscViewerDestroy(&vheader);
60: return(0);
61: }
63: static PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems)
64: {
65: PetscInt c, p;
66: GmshElement *elements;
67: int i, cellType, dim, numNodes, numNodesIgnore, numElem, numTags;
68: PetscInt pibuf[64];
69: int ibuf[16];
73: PetscMalloc1(numCells, &elements);
74: for (c = 0; c < numCells;) {
75: PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);
76: if (byteSwap) PetscByteSwap(&ibuf, PETSC_ENUM, 3);
77: if (binary) {
78: cellType = ibuf[0];
79: numElem = ibuf[1];
80: numTags = ibuf[2];
81: } else {
82: elements[c].id = ibuf[0];
83: cellType = ibuf[1];
84: numTags = ibuf[2];
85: numElem = 1;
86: }
87: /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */
88: numNodesIgnore = 0;
89: switch (cellType) {
90: case 1: /* 2-node line */
91: dim = 1;
92: numNodes = 2;
93: break;
94: case 2: /* 3-node triangle */
95: dim = 2;
96: numNodes = 3;
97: break;
98: case 3: /* 4-node quadrangle */
99: dim = 2;
100: numNodes = 4;
101: break;
102: case 4: /* 4-node tetrahedron */
103: dim = 3;
104: numNodes = 4;
105: break;
106: case 5: /* 8-node hexahedron */
107: dim = 3;
108: numNodes = 8;
109: break;
110: case 8: /* 3-node 2nd order line */
111: dim = 1;
112: numNodes = 2;
113: numNodesIgnore = 1;
114: break;
115: case 9: /* 6-node 2nd order triangle */
116: dim = 2;
117: numNodes = 3;
118: numNodesIgnore = 3;
119: break;
120: case 15: /* 1-node vertex */
121: dim = 0;
122: numNodes = 1;
123: break;
124: case 6: /* 6-node prism */
125: case 7: /* 5-node pyramid */
126: case 10: /* 9-node 2nd order quadrangle */
127: case 11: /* 10-node 2nd order tetrahedron */
128: case 12: /* 27-node 2nd order hexhedron */
129: case 13: /* 19-node 2nd order prism */
130: case 14: /* 14-node 2nd order pyramid */
131: default:
132: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
133: }
134: if (binary) {
135: const PetscInt nint = numNodes + numTags + 1 + numNodesIgnore;
136: for (i = 0; i < numElem; ++i, ++c) {
137: /* Loop over inner binary element block */
138: elements[c].dim = dim;
139: elements[c].numNodes = numNodes;
140: elements[c].numTags = numTags;
142: PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);
143: if (byteSwap) PetscByteSwap( &ibuf, PETSC_ENUM, nint);
144: elements[c].id = ibuf[0];
145: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
146: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
147: }
148: } else {
149: elements[c].dim = dim;
150: elements[c].numNodes = numNodes;
151: elements[c].numTags = numTags;
152: PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
153: PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
154: PetscViewerRead(viewer, pibuf, numNodesIgnore, NULL, PETSC_ENUM);
155: c++;
156: }
157: }
158: *gmsh_elems = elements;
159: return(0);
160: }
162: /*@
163: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
165: Collective on comm
167: Input Parameters:
168: + comm - The MPI communicator
169: . viewer - The Viewer associated with a Gmsh file
170: - interpolate - Create faces and edges in the mesh
172: Output Parameter:
173: . dm - The DM object representing the mesh
175: Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format
176: and http://www.geuz.org/gmsh/doc/texinfo/#MSH-binary-file-format
178: Level: beginner
180: .keywords: mesh,Gmsh
181: .seealso: DMPLEX, DMCreate()
182: @*/
183: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
184: {
185: PetscViewerType vtype;
186: GmshElement *gmsh_elem = NULL;
187: PetscSection coordSection;
188: Vec coordinates;
189: PetscBT periodicV = NULL, periodicC = NULL;
190: PetscScalar *coords;
191: PetscReal *coordsIn = NULL;
192: PetscInt dim = 0, embedDim, coordSize, c, v, d, r, cell, *periodicMap = NULL, *periodicMapI = NULL;
193: int i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum, shift = 1;
194: PetscMPIInt num_proc, rank;
195: char line[PETSC_MAX_PATH_LEN];
196: PetscBool zerobase = PETSC_FALSE, isbd = PETSC_FALSE, match, binary, bswap = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE;
200: PetscOptionsGetBool(NULL, NULL, "-gmsh_zero_base", &zerobase, NULL);
201: if (zerobase) shift = 0;
202: MPI_Comm_rank(comm, &rank);
203: MPI_Comm_size(comm, &num_proc);
204: DMCreate(comm, dm);
205: DMSetType(*dm, DMPLEX);
206: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
207: PetscViewerGetType(viewer, &vtype);
208: PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);
209: PetscOptionsGetBool(NULL,NULL,"-dm_plex_gmsh_periodic",&periodic,NULL);
210: PetscOptionsGetBool(NULL,NULL,"-dm_plex_gmsh_usemarker",&usemarker,NULL);
211: if (!rank || binary) {
212: PetscBool match;
213: int fileType, dataSize;
214: float version;
216: /* Read header */
217: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
218: PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
219: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
220: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
221: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
222: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
223: if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
224: if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
225: if (binary) {
226: int checkInt;
227: PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
228: if (checkInt != 1) {
229: PetscByteSwap(&checkInt, PETSC_ENUM, 1);
230: if (checkInt == 1) bswap = PETSC_TRUE;
231: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
232: }
233: } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
234: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
235: PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
236: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
237: /* OPTIONAL Read physical names */
238: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
239: PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);
240: if (match) {
241: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
242: snum = sscanf(line, "%d", &numRegions);
243: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
244: for (r = 0; r < numRegions; ++r) {PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);}
245: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
246: PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);
247: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
248: /* Initial read for vertex section */
249: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
250: }
251: /* Read vertices */
252: PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
253: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
254: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
255: snum = sscanf(line, "%d", &numVertices);
256: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
257: PetscMalloc1(numVertices*3, &coordsIn);
258: if (binary) {
259: size_t doubleSize, intSize;
260: PetscInt elementSize;
261: char *buffer;
262: double *baseptr;
263: PetscDataTypeGetSize(PETSC_ENUM, &intSize);
264: PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);
265: elementSize = (intSize + 3*doubleSize);
266: PetscMalloc1(elementSize*numVertices, &buffer);
267: PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);
268: if (bswap) PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);
269: for (v = 0; v < numVertices; ++v) {
270: baseptr = ((double*)(buffer+v*elementSize+intSize));
271: coordsIn[v*3+0] = (PetscReal) baseptr[0];
272: coordsIn[v*3+1] = (PetscReal) baseptr[1];
273: coordsIn[v*3+2] = (PetscReal) baseptr[2];
274: }
275: PetscFree(buffer);
276: } else {
277: for (v = 0; v < numVertices; ++v) {
278: PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);
279: PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);
280: if (i != (int) v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+shift);
281: }
282: }
283: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
284: PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
285: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
286: /* Read cells */
287: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
288: PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
289: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
290: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
291: snum = sscanf(line, "%d", &numCells);
292: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
293: }
295: if (!rank || binary) {
296: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
297: file contents multiple times to figure out the true number of cells and facets
298: in the given mesh. To make this more efficient we read the file contents only
299: once and store them in memory, while determining the true number of cells. */
300: DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);
301: for (trueNumCells=0, c = 0; c < numCells; ++c) {
302: if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
303: if (gmsh_elem[c].dim == dim) trueNumCells++;
304: }
305: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
306: PetscStrncmp(line, "$EndElements", 12, &match);
307: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
308: if (periodic) {
309: PetscInt pVert, *periodicMapT, *aux;
310: int numPeriodic;
312: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
313: PetscStrncmp(line, "$Periodic", 9, &match);
314: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
315: PetscMalloc1(numVertices, &periodicMapT);
316: PetscBTCreate(numVertices, &periodicV);
317: for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
318: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
319: snum = sscanf(line, "%d", &numPeriodic);
320: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
321: for (i = 0; i < numPeriodic; i++) {
322: int j, edim = -1, slaveTag = -1, masterTag = -1, nNodes;
324: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
325: snum = sscanf(line, "%d %d %d", &edim, &slaveTag, &masterTag); /* slaveTag and masterTag are unused */
326: if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
327: /* Type of tranformation, discarded */
328: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
329: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
330: snum = sscanf(line, "%d", &nNodes);
331: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
332: for (j = 0; j < nNodes; j++) {
333: PetscInt slaveNode, masterNode;
335: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
336: snum = sscanf(line, "%" PetscInt_FMT " %" PetscInt_FMT, &slaveNode, &masterNode);
337: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
338: periodicMapT[slaveNode - shift] = masterNode - shift;
339: PetscBTSet(periodicV, slaveNode - shift);
340: PetscBTSet(periodicV, masterNode - shift);
341: }
342: }
343: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
344: PetscStrncmp(line, "$EndPeriodic", 12, &match);
345: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
346: /* we may have slaves of slaves */
347: for (i = 0; i < numVertices; i++) {
348: while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
349: periodicMapT[i] = periodicMapT[periodicMapT[i]];
350: }
351: }
352: /* periodicMap : from old to new numbering (periodic vertices excluded)
353: periodicMapI: from new to old numbering */
354: PetscMalloc1(numVertices, &periodicMap);
355: PetscMalloc1(numVertices, &periodicMapI);
356: PetscMalloc1(numVertices, &aux);
357: for (i = 0, pVert = 0; i < numVertices; i++) {
358: if (periodicMapT[i] != i) {
359: pVert++;
360: } else {
361: aux[i] = i - pVert;
362: periodicMapI[i - pVert] = i;
363: }
364: }
365: for (i = 0 ; i < numVertices; i++) {
366: periodicMap[i] = aux[periodicMapT[i]];
367: }
368: PetscFree(periodicMapT);
369: PetscFree(aux);
370: /* remove periodic vertices */
371: numVertices = numVertices - pVert;
372: }
373: }
374: /* For binary we read on all ranks, but only build the plex on rank 0 */
375: if (binary && rank) {trueNumCells = 0; numVertices = 0;};
376: /* Allocate the cell-vertex mesh */
377: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
378: if (!rank) {
379: for (cell = 0, c = 0; c < numCells; ++c) {
380: if (gmsh_elem[c].dim == dim) {
381: DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);
382: cell++;
383: }
384: }
385: }
386: DMSetUp(*dm);
387: /* Add cell-vertex connections */
388: if (!rank) {
389: PetscInt pcone[8], corner;
391: for (cell = 0, c = 0; c < numCells; ++c) {
392: if (gmsh_elem[c].dim == dim) {
393: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
394: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
395: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
396: }
397: if (dim == 3) {
398: /* Tetrahedra are inverted */
399: if (gmsh_elem[c].numNodes == 4) {
400: PetscInt tmp = pcone[0];
401: pcone[0] = pcone[1];
402: pcone[1] = tmp;
403: }
404: /* Hexahedra are inverted */
405: if (gmsh_elem[c].numNodes == 8) {
406: PetscInt tmp = pcone[1];
407: pcone[1] = pcone[3];
408: pcone[3] = tmp;
409: }
410: }
411: DMPlexSetCone(*dm, cell, pcone);
412: cell++;
413: }
414: }
415: }
416: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
417: DMSetDimension(*dm, dim);
418: DMPlexSymmetrize(*dm);
419: DMPlexStratify(*dm);
420: if (interpolate) {
421: DM idm = NULL;
423: DMPlexInterpolate(*dm, &idm);
424: DMDestroy(dm);
425: *dm = idm;
426: }
428: if (!rank && usemarker) {
429: PetscInt f, fStart, fEnd;
431: DMCreateLabel(*dm, "marker");
432: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
433: for (f = fStart; f < fEnd; ++f) {
434: PetscInt suppSize;
436: DMPlexGetSupportSize(*dm, f, &suppSize);
437: if (suppSize == 1) {
438: PetscInt *cone = NULL, coneSize, p;
440: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
441: for (p = 0; p < coneSize; p += 2) {
442: DMSetLabelValue(*dm, "marker", cone[p], 1);
443: }
444: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
445: }
446: }
447: }
449: if (!rank) {
450: /* Apply boundary IDs by finding the relevant facets with vertex joins */
451: PetscInt pcone[8], corner, vStart, vEnd;
453: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
454: for (c = 0; c < numCells; ++c) {
455: if (gmsh_elem[c].dim == dim-1) {
456: const PetscInt *join;
457: PetscInt joinSize;
459: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
460: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
461: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
462: }
463: DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
464: if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
465: DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
466: DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
467: }
468: }
470: /* Create cell sets */
471: for (cell = 0, c = 0; c < numCells; ++c) {
472: if (gmsh_elem[c].dim == dim) {
473: if (gmsh_elem[c].numTags > 0) {
474: DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);
475: }
476: cell++;
477: }
478: }
480: /* Create vertex sets */
481: for (c = 0; c < numCells; ++c) {
482: if (gmsh_elem[c].dim == 0) {
483: if (gmsh_elem[c].numTags > 0) {
484: const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
485: PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
486: DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
487: }
488: }
489: }
490: }
492: /* Read coordinates */
493: embedDim = dim;
494: PetscOptionsGetBool(NULL, NULL, "-gmsh_bd", &isbd, NULL);
495: if (isbd) embedDim = dim+1;
496: DMGetCoordinateSection(*dm, &coordSection);
497: PetscSectionSetNumFields(coordSection, 1);
498: PetscSectionSetFieldComponents(coordSection, 0, embedDim);
499: if (periodicMap) { /* we need to localize coordinates on cells */
500: PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
501: } else {
502: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
503: }
504: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
505: PetscSectionSetDof(coordSection, v, embedDim);
506: PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
507: }
508: if (!rank && periodicMap) {
509: PetscBTCreate(trueNumCells, &periodicC);
510: for (cell = 0, c = 0; c < numCells; ++c) {
511: if (gmsh_elem[c].dim == dim) {
512: PetscBool pc = PETSC_FALSE;
513: PetscInt corner;
514: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
515: pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
516: }
517: if (pc) {
518: PetscBTSet(periodicC, cell);
519: PetscSectionSetDof(coordSection, cell, gmsh_elem[c].numNodes*dim);
520: PetscSectionSetFieldDof(coordSection, cell, 0, gmsh_elem[c].numNodes*dim);
521: }
522: cell++;
523: }
524: }
525: }
526: PetscSectionSetUp(coordSection);
527: PetscSectionGetStorageSize(coordSection, &coordSize);
528: VecCreate(PETSC_COMM_SELF, &coordinates);
529: PetscObjectSetName((PetscObject) coordinates, "coordinates");
530: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
531: VecSetBlockSize(coordinates, embedDim);
532: VecSetType(coordinates, VECSTANDARD);
533: VecGetArray(coordinates, &coords);
534: if (!rank) {
536: if (periodicMap) {
537: PetscInt off;
539: for (cell = 0, c = 0; c < numCells; ++c) {
540: PetscInt pcone[8], corner;
541: if (gmsh_elem[c].dim == dim) {
542: if (PetscUnlikely(PetscBTLookup(periodicC,cell))) {
543: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
544: pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
545: }
546: if (dim == 3) {
547: /* Tetrahedra are inverted */
548: if (gmsh_elem[c].numNodes == 4) {
549: PetscInt tmp = pcone[0];
550: pcone[0] = pcone[1];
551: pcone[1] = tmp;
552: }
553: /* Hexahedra are inverted */
554: if (gmsh_elem[c].numNodes == 8) {
555: PetscInt tmp = pcone[1];
556: pcone[1] = pcone[3];
557: pcone[3] = tmp;
558: }
559: }
560: PetscSectionGetOffset(coordSection, cell, &off);
561: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
562: PetscInt v = pcone[corner];
563: for (d = 0; d < embedDim; ++d) {
564: coords[off++] = coordsIn[v*3+d];
565: }
566: }
567: }
568: cell++;
569: }
570: }
571: for (v = 0; v < numVertices; ++v) {
572: PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
573: for (d = 0; d < embedDim; ++d) {
574: coords[off+d] = coordsIn[periodicMapI[v]*3+d];
575: }
576: }
577: } else {
578: for (v = 0; v < numVertices; ++v) {
579: for (d = 0; d < embedDim; ++d) {
580: coords[v*embedDim+d] = coordsIn[v*3+d];
581: }
582: }
583: }
584: }
585: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
586: VecRestoreArray(coordinates, &coords);
587: PetscFree(coordsIn);
588: DMSetCoordinatesLocal(*dm, coordinates);
589: VecDestroy(&coordinates);
590: PetscFree(periodicMap);
591: PetscFree(periodicMapI);
592: PetscBTDestroy(&periodicV);
593: PetscBTDestroy(&periodicC);
594: /* Clean up intermediate storage */
595: if (!rank || binary) PetscFree(gmsh_elem);
596: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
597: return(0);
598: }