Actual source code: plexply.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>

  4: /*@C
  5:   DMPlexCreatePLYFromFile - Create a DMPlex mesh from a PLY file.

  7: + comm        - The MPI communicator
  8: . filename    - Name of the .med file
  9: - interpolate - Create faces and edges in the mesh

 11:   Output Parameter:
 12: . dm  - The DM object representing the mesh

 14:   Note: https://en.wikipedia.org/wiki/PLY_(file_format)

 16:   Level: beginner

 18: .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateMedFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
 19: @*/
 20: PetscErrorCode DMPlexCreatePLYFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 21: {
 22:   PetscViewer  viewer;
 23:   Vec          coordinates;
 24:   PetscSection coordSection;
 25:   PetscScalar *coords;
 26:   char         line[PETSC_MAX_PATH_LEN], ntype[16], itype[16], name[1024], vtype[16];
 27:   PetscBool    match, byteSwap = PETSC_FALSE;
 28:   PetscInt     dim = 2, cdim = 3, Nvp = 0, coordSize, xi = -1, yi = -1, zi = -1, v, c, p;
 29:   PetscMPIInt  rank;
 30:   int          snum, Nv, Nc;

 32:   MPI_Comm_rank(comm, &rank);
 33:   PetscViewerCreate(comm, &viewer);
 34:   PetscViewerSetType(viewer, PETSCVIEWERBINARY);
 35:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 36:   PetscViewerFileSetName(viewer, filename);
 37:   if (rank == 0) {
 38:     PetscBool isAscii, isBinaryBig, isBinaryLittle;

 40:     /* Check for PLY file */
 41:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 42:     PetscStrncmp(line, "ply", PETSC_MAX_PATH_LEN, &match);
 44:     /* Check PLY format */
 45:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 46:     PetscStrncmp(line, "format", PETSC_MAX_PATH_LEN, &match);
 48:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 49:     PetscStrncmp(line, "ascii", PETSC_MAX_PATH_LEN, &isAscii);
 50:     PetscStrncmp(line, "binary_big_endian", PETSC_MAX_PATH_LEN, &isBinaryBig);
 51:     PetscStrncmp(line, "binary_little_endian", PETSC_MAX_PATH_LEN, &isBinaryLittle);
 53:     if (isBinaryLittle) byteSwap = PETSC_TRUE;
 54:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 55:     PetscStrncmp(line, "1.0", PETSC_MAX_PATH_LEN, &match);
 57:     /* Ignore comments */
 58:     match = PETSC_TRUE;
 59:     while (match) {
 60:       char buf = '\0';
 61:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 62:       PetscStrncmp(line, "comment", PETSC_MAX_PATH_LEN, &match);
 63:       if (match)
 64:         while (buf != '\n') PetscViewerRead(viewer, &buf, 1, NULL, PETSC_CHAR);
 65:     }
 66:     /* Read vertex information */
 67:     PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match);
 69:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 70:     PetscStrncmp(line, "vertex", PETSC_MAX_PATH_LEN, &match);
 72:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 73:     snum = sscanf(line, "%d", &Nv);
 75:     match = PETSC_TRUE;
 76:     while (match) {
 77:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 78:       PetscStrncmp(line, "property", PETSC_MAX_PATH_LEN, &match);
 79:       if (match) {
 80:         PetscBool matchB;

 82:         PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
 84:         snum = sscanf(line, "%s %s", ntype, name);
 86:         PetscStrncmp(ntype, "float32", 16, &matchB);
 87:         if (matchB) {
 88:           vtype[Nvp] = 'f';
 89:         } else {
 90:           PetscStrncmp(ntype, "int32", 16, &matchB);
 91:           if (matchB) {
 92:             vtype[Nvp] = 'd';
 93:           } else {
 94:             PetscStrncmp(ntype, "uint8", 16, &matchB);
 95:             if (matchB) {
 96:               vtype[Nvp] = 'c';
 97:             } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse type in PLY file header: %s", line);
 98:           }
 99:         }
100:         PetscStrncmp(name, "x", 16, &matchB);
101:         if (matchB) xi = Nvp;
102:         PetscStrncmp(name, "y", 16, &matchB);
103:         if (matchB) yi = Nvp;
104:         PetscStrncmp(name, "z", 16, &matchB);
105:         if (matchB) zi = Nvp;
106:         ++Nvp;
107:       }
108:     }
109:     /* Read cell information */
110:     PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match);
112:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
113:     PetscStrncmp(line, "face", PETSC_MAX_PATH_LEN, &match);
115:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
116:     snum = sscanf(line, "%d", &Nc);
118:     PetscViewerRead(viewer, line, 5, NULL, PETSC_STRING);
119:     snum = sscanf(line, "property list %s %s %s", ntype, itype, name);
121:     PetscStrncmp(ntype, "uint8", 1024, &match);
123:     PetscStrncmp(name, "vertex_indices", 1024, &match);
125:     /* Header should terminate */
126:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
127:     PetscStrncmp(line, "end_header", PETSC_MAX_PATH_LEN, &match);
129:   } else {
130:     Nc = Nv = 0;
131:   }
132:   DMCreate(comm, dm);
133:   DMSetType(*dm, DMPLEX);
134:   DMPlexSetChart(*dm, 0, Nc + Nv);
135:   DMSetDimension(*dm, dim);
136:   DMSetCoordinateDim(*dm, cdim);
137:   /* Read coordinates */
138:   DMGetCoordinateSection(*dm, &coordSection);
139:   PetscSectionSetNumFields(coordSection, 1);
140:   PetscSectionSetFieldComponents(coordSection, 0, cdim);
141:   PetscSectionSetChart(coordSection, Nc, Nc + Nv);
142:   for (v = Nc; v < Nc + Nv; ++v) {
143:     PetscSectionSetDof(coordSection, v, cdim);
144:     PetscSectionSetFieldDof(coordSection, v, 0, cdim);
145:   }
146:   PetscSectionSetUp(coordSection);
147:   PetscSectionGetStorageSize(coordSection, &coordSize);
148:   VecCreate(PETSC_COMM_SELF, &coordinates);
149:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
150:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
151:   VecSetBlockSize(coordinates, cdim);
152:   VecSetType(coordinates, VECSTANDARD);
153:   VecGetArray(coordinates, &coords);
154:   if (rank == 0) {
155:     float rbuf[1];
156:     int   ibuf[1];

158:     for (v = 0; v < Nv; ++v) {
159:       for (p = 0; p < Nvp; ++p) {
160:         if (vtype[p] == 'f') {
161:           PetscViewerRead(viewer, &rbuf, 1, NULL, PETSC_FLOAT);
162:           if (byteSwap) PetscByteSwap(&rbuf, PETSC_FLOAT, 1);
163:           if (p == xi) coords[v * cdim + 0] = rbuf[0];
164:           else if (p == yi) coords[v * cdim + 1] = rbuf[0];
165:           else if (p == zi) coords[v * cdim + 2] = rbuf[0];
166:         } else if (vtype[p] == 'd') {
167:           PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_INT);
168:           if (byteSwap) PetscByteSwap(&ibuf, PETSC_INT, 1);
169:         } else if (vtype[p] == 'c') {
170:           PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);
171:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid vertex property type in PLY file");
172:       }
173:     }
174:   }
175:   VecRestoreArray(coordinates, &coords);
176:   DMSetCoordinatesLocal(*dm, coordinates);
177:   VecDestroy(&coordinates);
178:   /* Read topology */
179:   if (rank == 0) {
180:     char     ibuf[1];
181:     PetscInt vbuf[16], corners;

183:     /* Assume same cells */
184:     PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);
185:     corners = ibuf[0];
186:     for (c = 0; c < Nc; ++c) DMPlexSetConeSize(*dm, c, corners);
187:     DMSetUp(*dm);
188:     for (c = 0; c < Nc; ++c) {
189:       if (c > 0) PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);
191:       PetscViewerRead(viewer, &vbuf, ibuf[0], NULL, PETSC_INT);
192:       if (byteSwap) PetscByteSwap(&vbuf, PETSC_INT, ibuf[0]);
193:       for (v = 0; v < ibuf[0]; ++v) vbuf[v] += Nc;
194:       DMPlexSetCone(*dm, c, vbuf);
195:     }
196:   }
197:   DMPlexSymmetrize(*dm);
198:   DMPlexStratify(*dm);
199:   PetscViewerDestroy(&viewer);
200:   if (interpolate) {
201:     DM idm;

203:     DMPlexInterpolate(*dm, &idm);
204:     DMDestroy(dm);
205:     *dm = idm;
206:   }
207:   return 0;
208: }