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: }