Actual source code: plexply.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  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;
 31:   PetscErrorCode  ierr;

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

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

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

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

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

206:     DMPlexInterpolate(*dm, &idm);
207:     DMDestroy(dm);
208:     *dm  = idm;
209:   }
210:   return(0);
211: }