Actual source code: plexgmsh.c
petsc-3.5.4 2015-05-23
1: #define PETSCDM_DLL
2: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
6: /*@
7: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file.
9: Collective on comm
11: Input Parameters:
12: + comm - The MPI communicator
13: . viewer - The Viewer associated with a Gmsh file
14: - interpolate - Create faces and edges in the mesh
16: Output Parameter:
17: . dm - The DM object representing the mesh
19: Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format
21: Level: beginner
23: .keywords: mesh,Gmsh
24: .seealso: DMPLEX, DMCreate()
25: @*/
26: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
27: {
28: FILE *fd;
29: PetscSection coordSection;
30: Vec coordinates;
31: PetscScalar *coords, *coordsIn = NULL;
32: PetscInt dim = 0, tdim = 0, coordSize, c, v, d, numCorners;
33: int numVertices = 0, numCells = 0, trueNumCells = 0, numFacets = 0, cone[8], tags[2], cellNum, snum;
34: long fpos = 0;
35: PetscMPIInt num_proc, rank;
36: char line[PETSC_MAX_PATH_LEN];
37: PetscBool match;
41: MPI_Comm_rank(comm, &rank);
42: MPI_Comm_size(comm, &num_proc);
43: DMCreate(comm, dm);
44: DMSetType(*dm, DMPLEX);
45: if (!rank) {
46: PetscBool match;
47: int fileType, dataSize;
49: PetscViewerASCIIGetPointer(viewer, &fd);
50: /* Read header */
51: fgets(line, PETSC_MAX_PATH_LEN, fd);
52: PetscStrncmp(line, "$MeshFormat\n", PETSC_MAX_PATH_LEN, &match);
53: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
54: snum = fscanf(fd, "2.2 %d %d\n", &fileType, &dataSize);CHKERRQ(snum != 2);
55: if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
56: if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
57: fgets(line, PETSC_MAX_PATH_LEN, fd);
58: PetscStrncmp(line, "$EndMeshFormat\n", PETSC_MAX_PATH_LEN, &match);
59: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
60: /* Read vertices */
61: fgets(line, PETSC_MAX_PATH_LEN, fd);
62: PetscStrncmp(line, "$Nodes\n", PETSC_MAX_PATH_LEN, &match);
63: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
64: snum = fscanf(fd, "%d\n", &numVertices);CHKERRQ(snum != 1);
65: PetscMalloc(numVertices*3 * sizeof(PetscScalar), &coordsIn);
66: for (v = 0; v < numVertices; ++v) {
67: double x, y, z;
68: int i;
70: snum = fscanf(fd, "%d %lg %lg %lg\n", &i, &x, &y, &z);CHKERRQ(snum != 4);
71: coordsIn[v*3+0] = x; coordsIn[v*3+1] = y; coordsIn[v*3+2] = z;
72: if (i != v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1);
73: }
74: fgets(line, PETSC_MAX_PATH_LEN, fd);
75: PetscStrncmp(line, "$EndNodes\n", PETSC_MAX_PATH_LEN, &match);
76: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
77: /* Read cells */
78: fgets(line, PETSC_MAX_PATH_LEN, fd);
79: PetscStrncmp(line, "$Elements\n", PETSC_MAX_PATH_LEN, &match);
80: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
81: snum = fscanf(fd, "%d\n", &numCells);CHKERRQ(snum != 1);
82: }
84: if (!rank) {
85: fpos = ftell(fd);
86: /* The Gmsh format disguises facets as elements, so we have to run through all "element" entries
87: to get the correct numCells and decide the topological dimension of the mesh */
88: trueNumCells = 0;
89: for (c = 0; c < numCells; ++c) {
90: DMPlexCreateGmsh_ReadElement(fd, &dim, &cellNum, &numCorners, cone, tags);
91: if (dim > tdim) {
92: tdim = dim;
93: trueNumCells = 0;
94: }
95: trueNumCells++;
96: }
97: }
98: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
99: numFacets = numCells - trueNumCells;
100: if (!rank) {
101: fseek(fd, fpos, SEEK_SET);
102: for (c = 0; c < numCells; ++c) {
103: DMPlexCreateGmsh_ReadElement(fd, &dim, &cellNum, &numCorners, cone, tags);
104: if (dim == tdim) {
105: DMPlexSetConeSize(*dm, c-numFacets, numCorners);
106: }
107: if (cellNum != c+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell number %d should be %d", cellNum, c+1);
108: }
109: }
110: DMSetUp(*dm);
111: if (!rank) {
112: PetscInt pcone[8], corner;
114: fseek(fd, fpos, SEEK_SET);
115: for (c = 0; c < numCells; ++c) {
116: DMPlexCreateGmsh_ReadElement(fd, &dim, &cellNum, &numCorners, cone, tags);
117: if (dim == tdim) {
118: for (corner = 0; corner < numCorners; ++corner) pcone[corner] = cone[corner] + trueNumCells-1;
119: DMPlexSetCone(*dm, c-numFacets, (const PetscInt *) pcone);
120: }
121: if (cellNum != c+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell number %d should be %d", cellNum, c+1);
122: }
123: fgets(line, PETSC_MAX_PATH_LEN, fd);
124: PetscStrncmp(line, "$EndElements\n", PETSC_MAX_PATH_LEN, &match);
125: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
126: }
127: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
128: DMPlexSetDimension(*dm, dim);
129: DMPlexSymmetrize(*dm);
130: DMPlexStratify(*dm);
131: if (interpolate) {
132: DM idm = NULL;
134: DMPlexInterpolate(*dm, &idm);
135: DMDestroy(dm);
136: *dm = idm;
137: }
139: if (!rank) {
140: /* Apply boundary IDs by finding the relevant facets with vertex joins */
141: PetscInt pcone[8], corner, vStart, vEnd;
143: fseek(fd, fpos, SEEK_SET);
144: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
145: for (c = 0; c < numCells; ++c) {
146: DMPlexCreateGmsh_ReadElement(fd, &dim, &cellNum, &numCorners, cone, tags);
147: if (dim == tdim-1) {
148: PetscInt joinSize;
149: const PetscInt *join;
150: for (corner = 0; corner < numCorners; ++corner) pcone[corner] = cone[corner] + vStart - 1;
151: DMPlexGetFullJoin(*dm, numCorners, (const PetscInt *) pcone, &joinSize, &join);
152: if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", cellNum);
153: DMPlexSetLabelValue(*dm, "Face Sets", join[0], tags[0]);
154: DMPlexRestoreJoin(*dm, numCorners, (const PetscInt *) pcone, &joinSize, &join);
155: }
156: if (cellNum != c+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell number %d should be %d", cellNum, c+1);
157: }
158: fgets(line, PETSC_MAX_PATH_LEN, fd);
159: PetscStrncmp(line, "$EndElements\n", PETSC_MAX_PATH_LEN, &match);
160: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
161: }
163: /* Read coordinates */
164: DMGetCoordinateSection(*dm, &coordSection);
165: PetscSectionSetNumFields(coordSection, 1);
166: PetscSectionSetFieldComponents(coordSection, 0, dim);
167: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
168: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
169: PetscSectionSetDof(coordSection, v, dim);
170: PetscSectionSetFieldDof(coordSection, v, 0, dim);
171: }
172: PetscSectionSetUp(coordSection);
173: PetscSectionGetStorageSize(coordSection, &coordSize);
174: VecCreate(comm, &coordinates);
175: PetscObjectSetName((PetscObject) coordinates, "coordinates");
176: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
177: VecSetType(coordinates, VECSTANDARD);
178: VecGetArray(coordinates, &coords);
179: if (!rank) {
180: for (v = 0; v < numVertices; ++v) {
181: for (d = 0; d < dim; ++d) {
182: coords[v*dim+d] = coordsIn[v*3+d];
183: }
184: }
185: }
186: VecRestoreArray(coordinates, &coords);
187: PetscFree(coordsIn);
188: DMSetCoordinatesLocal(*dm, coordinates);
189: VecDestroy(&coordinates);
190: return(0);
191: }
195: PetscErrorCode DMPlexCreateGmsh_ReadElement(FILE *fd, PetscInt *dim, int *cellNum, PetscInt *numCorners, int cone[], int tags[])
196: {
197: PetscInt t;
198: int numTags, snum, cellType;
201: snum = fscanf(fd, "%d %d %d", cellNum, &cellType, &numTags);CHKERRQ(snum != 3);
202: for (t = 0; t < numTags; ++t) {snum = fscanf(fd, "%d", &tags[t]);CHKERRQ(snum != 1);}
203: switch (cellType) {
204: case 1: /* 2-node line */
205: *dim = 1;
206: *numCorners = 2;
207: snum = fscanf(fd, "%d %d\n", &cone[0], &cone[1]);CHKERRQ(snum != *numCorners);
208: break;
209: case 2: /* 3-node triangle */
210: *dim = 2;
211: *numCorners = 3;
212: snum = fscanf(fd, "%d %d %d\n", &cone[0], &cone[1], &cone[2]);CHKERRQ(snum != *numCorners);
213: break;
214: case 3: /* 4-node quadrangle */
215: *dim = 2;
216: *numCorners = 4;
217: snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != *numCorners);
218: break;
219: case 4: /* 4-node tetrahedron */
220: *dim = 3;
221: *numCorners = 4;
222: snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != *numCorners);
223: break;
224: case 5: /* 8-node hexahedron */
225: *dim = 3;
226: *numCorners = 8;
227: snum = fscanf(fd, "%d %d %d %d %d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3], &cone[4], &cone[5], &cone[6], &cone[7]);CHKERRQ(snum != *numCorners);
228: break;
229: default:
230: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
231: }
232: return(0);
233: }