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:     else 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) while (buf != '\n') PetscViewerRead(viewer, &buf, 1, NULL, PETSC_CHAR);
 64:     }
 65:     /* Read vertex information */
 66:     PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match);
 68:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 69:     PetscStrncmp(line, "vertex", PETSC_MAX_PATH_LEN, &match);
 71:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 72:     snum = sscanf(line, "%d", &Nv);
 74:     match = PETSC_TRUE;
 75:     while (match) {
 76:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
 77:       PetscStrncmp(line, "property", PETSC_MAX_PATH_LEN, &match);
 78:       if (match) {
 79:         PetscBool matchB;

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

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

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

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