Actual source code: plexgmsh.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  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: }