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:   Input Parameters:
  8: + comm        - The MPI communicator
  9: . filename    - Name of the .med file
 10: - interpolate - Create faces and edges in the mesh

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

 15:   References:
 16: . * - https://en.wikipedia.org/wiki/PLY_(file_format)

 18:   Level: beginner

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

 34:   PetscFunctionBegin;
 35:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 36:   PetscCall(PetscViewerCreate(comm, &viewer));
 37:   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
 38:   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
 39:   PetscCall(PetscViewerFileSetName(viewer, filename));
 40:   if (rank == 0) {
 41:     PetscBool isAscii, isBinaryBig, isBinaryLittle;

 43:     /* Check for PLY file */
 44:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 45:     PetscCall(PetscStrncmp(line, "ply", PETSC_MAX_PATH_LEN, &match));
 46:     PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file");
 47:     /* Check PLY format */
 48:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 49:     PetscCall(PetscStrncmp(line, "format", PETSC_MAX_PATH_LEN, &match));
 50:     PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file");
 51:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 52:     PetscCall(PetscStrncmp(line, "ascii", PETSC_MAX_PATH_LEN, &isAscii));
 53:     PetscCall(PetscStrncmp(line, "binary_big_endian", PETSC_MAX_PATH_LEN, &isBinaryBig));
 54:     PetscCall(PetscStrncmp(line, "binary_little_endian", PETSC_MAX_PATH_LEN, &isBinaryLittle));
 55:     PetscCheck(!isAscii, PETSC_COMM_SELF, PETSC_ERR_SUP, "PLY ascii format not yet supported");
 56:     if (isBinaryLittle) byteSwap = PETSC_TRUE;
 57:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 58:     PetscCall(PetscStrncmp(line, "1.0", PETSC_MAX_PATH_LEN, &match));
 59:     PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid version of PLY file, %s", line);
 60:     /* Ignore comments */
 61:     match = PETSC_TRUE;
 62:     while (match) {
 63:       char buf = '\0';
 64:       PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 65:       PetscCall(PetscStrncmp(line, "comment", PETSC_MAX_PATH_LEN, &match));
 66:       if (match)
 67:         while (buf != '\n') PetscCall(PetscViewerRead(viewer, &buf, 1, NULL, PETSC_CHAR));
 68:     }
 69:     /* Read vertex information */
 70:     PetscCall(PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match));
 71:     PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
 72:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 73:     PetscCall(PetscStrncmp(line, "vertex", PETSC_MAX_PATH_LEN, &match));
 74:     PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
 75:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 76:     snum = sscanf(line, "%d", &Nv);
 77:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
 78:     match = PETSC_TRUE;
 79:     while (match) {
 80:       PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
 81:       PetscCall(PetscStrncmp(line, "property", PETSC_MAX_PATH_LEN, &match));
 82:       if (match) {
 83:         PetscBool matchB;

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

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

186:     /* Assume same cells */
187:     PetscCall(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR));
188:     corners = ibuf[0];
189:     for (c = 0; c < Nc; ++c) PetscCall(DMPlexSetConeSize(*dm, c, corners));
190:     PetscCall(DMSetUp(*dm));
191:     for (c = 0; c < Nc; ++c) {
192:       if (c > 0) PetscCall(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR));
193:       PetscCheck(ibuf[0] == corners, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "All cells must have the same number of vertices in PLY file: %" PetscInt_FMT " != %" PetscInt_FMT, (PetscInt)ibuf[0], corners);
194:       PetscCall(PetscViewerRead(viewer, &vbuf, ibuf[0], NULL, PETSC_INT));
195:       if (byteSwap) PetscCall(PetscByteSwap(&vbuf, PETSC_INT, ibuf[0]));
196:       for (v = 0; v < ibuf[0]; ++v) vbuf[v] += Nc;
197:       PetscCall(DMPlexSetCone(*dm, c, vbuf));
198:     }
199:   }
200:   PetscCall(DMPlexSymmetrize(*dm));
201:   PetscCall(DMPlexStratify(*dm));
202:   PetscCall(PetscViewerDestroy(&viewer));
203:   if (interpolate) {
204:     DM idm;

206:     PetscCall(DMPlexInterpolate(*dm, &idm));
207:     PetscCall(DMDestroy(dm));
208:     *dm = idm;
209:   }
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }