Actual source code: plexply.c
petsc-3.9.4 2018-09-11
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: }