Actual source code: plexgmsh.c
petsc-3.9.4 2018-09-11
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 8: /* 3-node 2nd order line */
128: dim = 1;
129: numNodes = 2;
130: numNodesIgnore = 1;
131: break;
132: case 9: /* 6-node 2nd order triangle */
133: dim = 2;
134: numNodes = 3;
135: numNodesIgnore = 3;
136: break;
137: case 15: /* 1-node vertex */
138: dim = 0;
139: numNodes = 1;
140: break;
141: case 6: /* 6-node prism */
142: case 7: /* 5-node pyramid */
143: case 10: /* 9-node 2nd order quadrangle */
144: case 11: /* 10-node 2nd order tetrahedron */
145: case 12: /* 27-node 2nd order hexhedron */
146: case 13: /* 19-node 2nd order prism */
147: case 14: /* 14-node 2nd order pyramid */
148: default:
149: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
150: }
151: if (binary) {
152: const int nint = 1 + numTags + numNodes + numNodesIgnore;
153: /* Loop over element blocks */
154: for (i = 0; i < numElem; ++i, ++c) {
155: PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
156: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, nint);}
157: elements[c].dim = dim;
158: elements[c].numNodes = numNodes;
159: elements[c].numTags = numTags;
160: elements[c].id = ibuf[0];
161: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
162: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
163: }
164: } else {
165: elements[c].dim = dim;
166: elements[c].numNodes = numNodes;
167: elements[c].numTags = numTags;
168: PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
169: PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
170: PetscViewerRead(viewer, ibuf, numNodesIgnore, NULL, PETSC_ENUM);
171: c++;
172: }
173: }
174: *gmsh_elems = elements;
175: return(0);
176: }
178: /*@
179: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
181: Collective on comm
183: Input Parameters:
184: + comm - The MPI communicator
185: . viewer - The Viewer associated with a Gmsh file
186: - interpolate - Create faces and edges in the mesh
188: Output Parameter:
189: . dm - The DM object representing the mesh
191: Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format
192: and http://www.geuz.org/gmsh/doc/texinfo/#MSH-binary-file-format
194: Level: beginner
196: .keywords: mesh,Gmsh
197: .seealso: DMPLEX, DMCreate()
198: @*/
199: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
200: {
201: PetscViewer parentviewer = NULL;
202: double *coordsIn = NULL;
203: GmshElement *gmsh_elem = NULL;
204: PetscSection coordSection;
205: Vec coordinates;
206: PetscBT periodicV = NULL, periodicC = NULL;
207: PetscScalar *coords;
208: PetscInt dim = 0, embedDim, coordSize, c, v, d, r, cell, *periodicMap = NULL, *periodicMapI = NULL;
209: int i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum, shift = 1;
210: PetscMPIInt rank;
211: char line[PETSC_MAX_PATH_LEN];
212: PetscBool binary, byteSwap = PETSC_FALSE, zerobase = PETSC_FALSE, isbd = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE;
216: MPI_Comm_rank(comm, &rank);
217: PetscOptionsGetBool(NULL, NULL, "-dm_plex_gmsh_periodic", &periodic, NULL);
218: PetscOptionsGetBool(NULL, NULL, "-dm_plex_gmsh_use_marker", &usemarker, NULL);
219: PetscOptionsGetBool(NULL, NULL, "-gmsh_zero_base", &zerobase, NULL);
220: PetscOptionsGetBool(NULL, NULL, "-gmsh_bd", &isbd, NULL);
221: if (zerobase) shift = 0;
223: DMCreate(comm, dm);
224: DMSetType(*dm, DMPLEX);
225: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
227: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
229: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
230: if (binary) {
231: parentviewer = viewer;
232: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
233: }
235: if (!rank) {
236: PetscBool match;
237: int fileType, dataSize;
238: float version;
240: /* Read header */
241: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
242: PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
243: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
244: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
245: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
246: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
247: if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
248: if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
249: if (binary) {
250: int checkInt;
251: PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
252: if (checkInt != 1) {
253: PetscByteSwap(&checkInt, PETSC_ENUM, 1);
254: if (checkInt == 1) byteSwap = PETSC_TRUE;
255: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
256: }
257: } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
258: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
259: PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
260: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
262: /* OPTIONAL Read physical names */
263: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
264: PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);
265: if (match) {
266: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
267: snum = sscanf(line, "%d", &numRegions);
268: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
269: for (r = 0; r < numRegions; ++r) {
270: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
271: }
272: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
273: PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);
274: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
275: /* Initial read for vertex section */
276: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
277: }
279: /* Read vertices */
280: PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
281: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
282: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
283: snum = sscanf(line, "%d", &numVertices);
284: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
285: DMPlexCreateGmsh_ReadNodes(viewer, binary, byteSwap, shift, numVertices, &coordsIn);
286: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
287: PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
288: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
290: /* Read cells */
291: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
292: PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
293: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
294: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
295: snum = sscanf(line, "%d", &numCells);
296: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
297: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
298: file contents multiple times to figure out the true number of cells and facets
299: in the given mesh. To make this more efficient we read the file contents only
300: once and store them in memory, while determining the true number of cells. */
301: DMPlexCreateGmsh_ReadElements(viewer, binary, byteSwap, shift, numCells, &gmsh_elem);
302: for (trueNumCells=0, c = 0; c < numCells; ++c) {
303: if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
304: if (gmsh_elem[c].dim == dim) trueNumCells++;
305: }
306: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
307: PetscStrncmp(line, "$EndElements", 12, &match);
308: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
310: /* OPTIONAL Read periodic section */
311: if (periodic) {
312: PetscInt pVert, *periodicMapT, *aux;
313: int numPeriodic;
315: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
316: PetscStrncmp(line, "$Periodic", 9, &match);
317: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
318: PetscMalloc1(numVertices, &periodicMapT);
319: PetscBTCreate(numVertices, &periodicV);
320: for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
321: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
322: snum = sscanf(line, "%d", &numPeriodic);
323: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
324: for (i = 0; i < numPeriodic; i++) {
325: int j, edim = -1, slaveTag = -1, masterTag = -1, nNodes, slaveNode, masterNode;
327: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
328: snum = sscanf(line, "%d %d %d", &edim, &slaveTag, &masterTag); /* slaveTag and masterTag are unused */
329: if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
330: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
331: snum = sscanf(line, "%d", &nNodes);
332: if (snum != 1) { /* discard tranformation and try again */
333: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
334: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
335: snum = sscanf(line, "%d", &nNodes);
336: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
337: }
338: for (j = 0; j < nNodes; j++) {
339: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
340: snum = sscanf(line, "%d %d", &slaveNode, &masterNode);
341: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
342: periodicMapT[slaveNode - shift] = masterNode - shift;
343: PetscBTSet(periodicV, slaveNode - shift);
344: PetscBTSet(periodicV, masterNode - shift);
345: }
346: }
347: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
348: PetscStrncmp(line, "$EndPeriodic", 12, &match);
349: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
350: /* we may have slaves of slaves */
351: for (i = 0; i < numVertices; i++) {
352: while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
353: periodicMapT[i] = periodicMapT[periodicMapT[i]];
354: }
355: }
356: /* periodicMap : from old to new numbering (periodic vertices excluded)
357: periodicMapI: from new to old numbering */
358: PetscMalloc1(numVertices, &periodicMap);
359: PetscMalloc1(numVertices, &periodicMapI);
360: PetscMalloc1(numVertices, &aux);
361: for (i = 0, pVert = 0; i < numVertices; i++) {
362: if (periodicMapT[i] != i) {
363: pVert++;
364: } else {
365: aux[i] = i - pVert;
366: periodicMapI[i - pVert] = i;
367: }
368: }
369: for (i = 0 ; i < numVertices; i++) {
370: periodicMap[i] = aux[periodicMapT[i]];
371: }
372: PetscFree(periodicMapT);
373: PetscFree(aux);
374: /* remove periodic vertices */
375: numVertices = numVertices - pVert;
376: }
377: }
379: if (parentviewer) {
380: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
381: }
383: /* Allocate the cell-vertex mesh */
384: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
385: for (cell = 0, c = 0; c < numCells; ++c) {
386: if (gmsh_elem[c].dim == dim) {
387: DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);
388: cell++;
389: }
390: }
391: DMSetUp(*dm);
392: /* Add cell-vertex connections */
393: for (cell = 0, c = 0; c < numCells; ++c) {
394: if (gmsh_elem[c].dim == dim) {
395: PetscInt pcone[8], corner;
396: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
397: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
398: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
399: }
400: if (dim == 3) {
401: /* Tetrahedra are inverted */
402: if (gmsh_elem[c].numNodes == 4) {
403: PetscInt tmp = pcone[0];
404: pcone[0] = pcone[1];
405: pcone[1] = tmp;
406: }
407: /* Hexahedra are inverted */
408: if (gmsh_elem[c].numNodes == 8) {
409: PetscInt tmp = pcone[1];
410: pcone[1] = pcone[3];
411: pcone[3] = tmp;
412: }
413: }
414: DMPlexSetCone(*dm, cell, pcone);
415: cell++;
416: }
417: }
418: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
419: DMSetDimension(*dm, dim);
420: DMPlexSymmetrize(*dm);
421: DMPlexStratify(*dm);
422: if (interpolate) {
423: DM idm;
425: DMPlexInterpolate(*dm, &idm);
426: DMDestroy(dm);
427: *dm = idm;
428: }
430: if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
431: if (!rank && usemarker) {
432: PetscInt f, fStart, fEnd;
434: DMCreateLabel(*dm, "marker");
435: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
436: for (f = fStart; f < fEnd; ++f) {
437: PetscInt suppSize;
439: DMPlexGetSupportSize(*dm, f, &suppSize);
440: if (suppSize == 1) {
441: PetscInt *cone = NULL, coneSize, p;
443: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
444: for (p = 0; p < coneSize; p += 2) {
445: DMSetLabelValue(*dm, "marker", cone[p], 1);
446: }
447: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
448: }
449: }
450: }
452: if (!rank) {
453: PetscInt vStart, vEnd;
455: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
456: for (cell = 0, c = 0; c < numCells; ++c) {
458: /* Create face sets */
459: if (interpolate && gmsh_elem[c].dim == dim-1) {
460: const PetscInt *join;
461: PetscInt joinSize, pcone[8], corner;
462: /* Find the relevant facet with vertex joins */
463: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
464: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
465: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
466: }
467: DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
468: if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
469: DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
470: DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
471: }
473: /* Create cell sets */
474: if (gmsh_elem[c].dim == dim) {
475: if (gmsh_elem[c].numTags > 0) {
476: DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);
477: }
478: cell++;
479: }
481: /* Create vertex sets */
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: const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
486: DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
487: }
488: }
489: }
490: }
492: /* Create coordinates */
493: embedDim = isbd ? dim+1 : dim;
494: DMGetCoordinateSection(*dm, &coordSection);
495: PetscSectionSetNumFields(coordSection, 1);
496: PetscSectionSetFieldComponents(coordSection, 0, embedDim);
497: if (periodicMap) { /* we need to localize coordinates on cells */
498: PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
499: } else {
500: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
501: }
502: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
503: PetscSectionSetDof(coordSection, v, embedDim);
504: PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
505: }
506: if (periodicMap) {
507: PetscBTCreate(trueNumCells, &periodicC);
508: for (cell = 0, c = 0; c < numCells; ++c) {
509: if (gmsh_elem[c].dim == dim) {
510: PetscInt corner;
511: PetscBool pc = PETSC_FALSE;
512: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
513: pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
514: }
515: if (pc) {
516: PetscInt dof = gmsh_elem[c].numNodes*embedDim;
517: PetscBTSet(periodicC, cell);
518: PetscSectionSetDof(coordSection, cell, dof);
519: PetscSectionSetFieldDof(coordSection, cell, 0, dof);
520: }
521: cell++;
522: }
523: }
524: }
525: PetscSectionSetUp(coordSection);
526: PetscSectionGetStorageSize(coordSection, &coordSize);
527: VecCreate(PETSC_COMM_SELF, &coordinates);
528: PetscObjectSetName((PetscObject) coordinates, "coordinates");
529: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
530: VecSetBlockSize(coordinates, embedDim);
531: VecSetType(coordinates, VECSTANDARD);
532: VecGetArray(coordinates, &coords);
533: if (periodicMap) {
534: PetscInt off;
536: for (cell = 0, c = 0; c < numCells; ++c) {
537: PetscInt pcone[8], corner;
538: if (gmsh_elem[c].dim == dim) {
539: if (PetscUnlikely(PetscBTLookup(periodicC, cell))) {
540: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
541: pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
542: }
543: if (dim == 3) {
544: /* Tetrahedra are inverted */
545: if (gmsh_elem[c].numNodes == 4) {
546: PetscInt tmp = pcone[0];
547: pcone[0] = pcone[1];
548: pcone[1] = tmp;
549: }
550: /* Hexahedra are inverted */
551: if (gmsh_elem[c].numNodes == 8) {
552: PetscInt tmp = pcone[1];
553: pcone[1] = pcone[3];
554: pcone[3] = tmp;
555: }
556: }
557: PetscSectionGetOffset(coordSection, cell, &off);
558: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
559: v = pcone[corner];
560: for (d = 0; d < embedDim; ++d) {
561: coords[off++] = (PetscReal) coordsIn[v*3+d];
562: }
563: }
564: }
565: cell++;
566: }
567: }
568: for (v = 0; v < numVertices; ++v) {
569: PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
570: for (d = 0; d < embedDim; ++d) {
571: coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
572: }
573: }
574: } else {
575: for (v = 0; v < numVertices; ++v) {
576: for (d = 0; d < embedDim; ++d) {
577: coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
578: }
579: }
580: }
581: VecRestoreArray(coordinates, &coords);
582: DMSetCoordinatesLocal(*dm, coordinates);
583: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
585: PetscFree(coordsIn);
586: PetscFree(gmsh_elem);
587: PetscFree(periodicMap);
588: PetscFree(periodicMapI);
589: PetscBTDestroy(&periodicV);
590: PetscBTDestroy(&periodicC);
591: VecDestroy(&coordinates);
593: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
594: return(0);
595: }