Actual source code: plexgmsh.c
petsc-3.6.1 2015-08-06
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
6: /*@C
7: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
9: + comm - The MPI communicator
10: . filename - Name of the Gmsh file
11: - interpolate - Create faces and edges in the mesh
13: Output Parameter:
14: . dm - The DM object representing the mesh
16: Level: beginner
18: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
19: @*/
20: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
21: {
22: PetscViewer viewer, vheader;
23: PetscMPIInt rank;
24: PetscViewerType vtype;
25: char line[PETSC_MAX_PATH_LEN];
26: int snum;
27: PetscBool match;
28: int fT;
29: PetscInt fileType;
30: float version;
31: PetscErrorCode ierr;
34: MPI_Comm_rank(comm, &rank);
35: /* Determine Gmsh file type (ASCII or binary) from file header */
36: PetscViewerCreate(comm, &vheader);
37: PetscViewerSetType(vheader, PETSCVIEWERASCII);
38: PetscViewerFileSetMode(vheader, FILE_MODE_READ);
39: PetscViewerFileSetName(vheader, filename);
40: if (!rank) {
41: /* Read only the first two lines of the Gmsh file */
42: PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);
43: PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &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, &fT);
47: fileType = (PetscInt) fT;
48: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
49: if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
50: }
51: MPI_Bcast(&fileType, 1, MPIU_INT, 0, comm);
52: /* Create appropriate viewer and build plex */
53: if (fileType == 0) vtype = PETSCVIEWERASCII;
54: else vtype = PETSCVIEWERBINARY;
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: PetscViewerDestroy(&vheader);
62: return(0);
63: }
68: /*@
69: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
71: Collective on comm
73: Input Parameters:
74: + comm - The MPI communicator
75: . viewer - The Viewer associated with a Gmsh file
76: - interpolate - Create faces and edges in the mesh
78: Output Parameter:
79: . dm - The DM object representing the mesh
81: Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format
82: and http://www.geuz.org/gmsh/doc/texinfo/#MSH-binary-file-format
84: Level: beginner
86: .keywords: mesh,Gmsh
87: .seealso: DMPLEX, DMCreate()
88: @*/
89: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
90: {
91: PetscViewerType vtype;
92: GmshElement *gmsh_elem;
93: PetscSection coordSection;
94: Vec coordinates;
95: PetscScalar *coords, *coordsIn = NULL;
96: PetscInt dim = 0, coordSize, c, v, d, cell;
97: int i, numVertices = 0, numCells = 0, trueNumCells = 0, snum;
98: PetscMPIInt num_proc, rank;
99: char line[PETSC_MAX_PATH_LEN];
100: PetscBool match, binary, bswap = PETSC_FALSE;
104: MPI_Comm_rank(comm, &rank);
105: MPI_Comm_size(comm, &num_proc);
106: DMCreate(comm, dm);
107: DMSetType(*dm, DMPLEX);
108: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
109: PetscViewerGetType(viewer, &vtype);
110: PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);
111: if (!rank || binary) {
112: PetscBool match;
113: int fileType, dataSize;
114: float version;
116: /* Read header */
117: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
118: PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
119: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
120: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
121: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
122: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
123: if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
124: if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
125: if (binary) {
126: int checkInt;
127: PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
128: if (checkInt != 1) {
129: PetscByteSwap(&checkInt, PETSC_ENUM, 1);
130: if (checkInt == 1) bswap = PETSC_TRUE;
131: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
132: }
133: } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
134: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
135: PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
136: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
137: /* Read vertices */
138: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
139: PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
140: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
141: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
142: snum = sscanf(line, "%d", &numVertices);
143: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
144: PetscMalloc1(numVertices*3, &coordsIn);
145: if (binary) {
146: size_t doubleSize, intSize;
147: PetscInt elementSize;
148: char *buffer;
149: PetscScalar *baseptr;
150: PetscDataTypeGetSize(PETSC_ENUM, &intSize);
151: PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);
152: elementSize = (intSize + 3*doubleSize);
153: PetscMalloc1(elementSize*numVertices, &buffer);
154: PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);
155: if (bswap) PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);
156: for (v = 0; v < numVertices; ++v) {
157: baseptr = ((PetscScalar*)(buffer+v*elementSize+intSize));
158: coordsIn[v*3+0] = baseptr[0];
159: coordsIn[v*3+1] = baseptr[1];
160: coordsIn[v*3+2] = baseptr[2];
161: }
162: PetscFree(buffer);
163: } else {
164: for (v = 0; v < numVertices; ++v) {
165: PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);
166: PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);
167: if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1);
168: }
169: }
170: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
171: PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
172: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
173: /* Read cells */
174: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
175: PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
176: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
177: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
178: snum = sscanf(line, "%d", &numCells);
179: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
180: }
182: if (!rank || binary) {
183: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
184: file contents multiple times to figure out the true number of cells and facets
185: in the given mesh. To make this more efficient we read the file contents only
186: once and store them in memory, while determining the true number of cells. */
187: DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);
188: for (trueNumCells=0, c = 0; c < numCells; ++c) {
189: if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
190: if (gmsh_elem[c].dim == dim) trueNumCells++;
191: }
192: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
193: PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);
194: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
195: }
196: /* For binary we read on all ranks, but only build the plex on rank 0 */
197: if (binary && rank) {trueNumCells = 0; numVertices = 0;};
198: /* Allocate the cell-vertex mesh */
199: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
200: if (!rank) {
201: for (cell = 0, c = 0; c < numCells; ++c) {
202: if (gmsh_elem[c].dim == dim) {
203: DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);
204: cell++;
205: }
206: }
207: }
208: DMSetUp(*dm);
209: /* Add cell-vertex connections */
210: if (!rank) {
211: PetscInt pcone[8], corner;
212: for (cell = 0, c = 0; c < numCells; ++c) {
213: if (gmsh_elem[c].dim == dim) {
214: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
215: pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1;
216: }
217: DMPlexSetCone(*dm, cell, pcone);
218: cell++;
219: }
220: }
221: }
222: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
223: DMSetDimension(*dm, dim);
224: DMPlexSymmetrize(*dm);
225: DMPlexStratify(*dm);
226: if (interpolate) {
227: DM idm = NULL;
229: DMPlexInterpolate(*dm, &idm);
230: DMDestroy(dm);
231: *dm = idm;
232: }
234: if (!rank) {
235: /* Apply boundary IDs by finding the relevant facets with vertex joins */
236: PetscInt pcone[8], corner, vStart, vEnd;
238: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
239: for (c = 0; c < numCells; ++c) {
240: if (gmsh_elem[c].dim == dim-1) {
241: PetscInt joinSize;
242: const PetscInt *join;
243: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
244: pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1;
245: }
246: DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
247: if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
248: DMPlexSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
249: DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
250: }
251: }
252: }
254: /* Read coordinates */
255: DMGetCoordinateSection(*dm, &coordSection);
256: PetscSectionSetNumFields(coordSection, 1);
257: PetscSectionSetFieldComponents(coordSection, 0, dim);
258: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
259: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
260: PetscSectionSetDof(coordSection, v, dim);
261: PetscSectionSetFieldDof(coordSection, v, 0, dim);
262: }
263: PetscSectionSetUp(coordSection);
264: PetscSectionGetStorageSize(coordSection, &coordSize);
265: VecCreate(comm, &coordinates);
266: PetscObjectSetName((PetscObject) coordinates, "coordinates");
267: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
268: VecSetType(coordinates, VECSTANDARD);
269: VecGetArray(coordinates, &coords);
270: if (!rank) {
271: for (v = 0; v < numVertices; ++v) {
272: for (d = 0; d < dim; ++d) {
273: coords[v*dim+d] = coordsIn[v*3+d];
274: }
275: }
276: }
277: VecRestoreArray(coordinates, &coords);
278: PetscFree(coordsIn);
279: DMSetCoordinatesLocal(*dm, coordinates);
280: VecDestroy(&coordinates);
281: /* Clean up intermediate storage */
282: if (!rank || binary) PetscFree(gmsh_elem);
283: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
284: return(0);
285: }
289: PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems)
290: {
291: PetscInt c, p;
292: GmshElement *elements;
293: int i, cellType, dim, numNodes, numElem, numTags;
294: int ibuf[16];
298: PetscMalloc1(numCells, &elements);
299: for (c = 0; c < numCells;) {
300: PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);
301: if (byteSwap) PetscByteSwap(&ibuf, PETSC_ENUM, 3);
302: if (binary) {
303: cellType = ibuf[0];
304: numElem = ibuf[1];
305: numTags = ibuf[2];
306: } else {
307: elements[c].id = ibuf[0];
308: cellType = ibuf[1];
309: numTags = ibuf[2];
310: numElem = 1;
311: }
312: switch (cellType) {
313: case 1: /* 2-node line */
314: dim = 1;
315: numNodes = 2;
316: break;
317: case 2: /* 3-node triangle */
318: dim = 2;
319: numNodes = 3;
320: break;
321: case 3: /* 4-node quadrangle */
322: dim = 2;
323: numNodes = 4;
324: break;
325: case 4: /* 4-node tetrahedron */
326: dim = 3;
327: numNodes = 4;
328: break;
329: case 5: /* 8-node hexahedron */
330: dim = 3;
331: numNodes = 8;
332: break;
333: case 15: /* 1-node vertex */
334: dim = 0;
335: numNodes = 1;
336: break;
337: default:
338: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
339: }
340: if (binary) {
341: const PetscInt nint = numNodes + numTags + 1;
342: for (i = 0; i < numElem; ++i, ++c) {
343: /* Loop over inner binary element block */
344: elements[c].dim = dim;
345: elements[c].numNodes = numNodes;
346: elements[c].numTags = numTags;
348: PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);
349: if (byteSwap) PetscByteSwap( &ibuf, PETSC_ENUM, nint);
350: elements[c].id = ibuf[0];
351: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
352: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
353: }
354: } else {
355: elements[c].dim = dim;
356: elements[c].numNodes = numNodes;
357: elements[c].numTags = numTags;
358: PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
359: PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
360: c++;
361: }
362: }
363: *gmsh_elems = elements;
364: return(0);
365: }