1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
3 /*@
4 DMPlexCreatePLYFromFile - Create a `DMPLEX` mesh from a PLY <https://en.wikipedia.org/wiki/PLY_(file_format)> file.
5
6 Input Parameters:
7 + comm - The MPI communicator
8 . filename - Name of the .ply file
9 - interpolate - Create faces and edges in the mesh
10
11 Output Parameter:
12 . dm - The `DMPLEX` object representing the mesh
13
14 Level: beginner
15
16 .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
17 @*/
DMPlexCreatePLYFromFile(MPI_Comm comm,const char filename[],PetscBool interpolate,DM * dm)18 PetscErrorCode DMPlexCreatePLYFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
19 {
20 PetscViewer viewer;
21 Vec coordinates;
22 PetscSection coordSection;
23 PetscScalar *coords;
24 char line[PETSC_MAX_PATH_LEN], ntype[16], itype[16], name[1024], vtype[16];
25 PetscBool match, byteSwap = PETSC_FALSE;
26 PetscInt dim = 2, cdim = 3, Nvp = 0, coordSize, xi = -1, yi = -1, zi = -1, v, c, p;
27 PetscMPIInt rank;
28 int snum, Nv, Nc;
29
30 PetscFunctionBegin;
31 PetscCallMPI(MPI_Comm_rank(comm, &rank));
32 PetscCall(PetscViewerCreate(comm, &viewer));
33 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
34 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
35 PetscCall(PetscViewerFileSetName(viewer, filename));
36 if (rank == 0) {
37 PetscBool isAscii, isBinaryBig, isBinaryLittle;
38
39 /* Check for PLY file */
40 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
41 PetscCall(PetscStrncmp(line, "ply", PETSC_MAX_PATH_LEN, &match));
42 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file");
43 /* Check PLY format */
44 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
45 PetscCall(PetscStrncmp(line, "format", PETSC_MAX_PATH_LEN, &match));
46 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file");
47 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
48 PetscCall(PetscStrncmp(line, "ascii", PETSC_MAX_PATH_LEN, &isAscii));
49 PetscCall(PetscStrncmp(line, "binary_big_endian", PETSC_MAX_PATH_LEN, &isBinaryBig));
50 PetscCall(PetscStrncmp(line, "binary_little_endian", PETSC_MAX_PATH_LEN, &isBinaryLittle));
51 PetscCheck(!isAscii, PETSC_COMM_SELF, PETSC_ERR_SUP, "PLY ascii format not yet supported");
52 if (isBinaryLittle) byteSwap = PETSC_TRUE;
53 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
54 PetscCall(PetscStrncmp(line, "1.0", PETSC_MAX_PATH_LEN, &match));
55 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid version of PLY file, %s", line);
56 /* Ignore comments */
57 match = PETSC_TRUE;
58 while (match) {
59 char buf = '\0';
60 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
61 PetscCall(PetscStrncmp(line, "comment", PETSC_MAX_PATH_LEN, &match));
62 if (match)
63 while (buf != '\n') PetscCall(PetscViewerRead(viewer, &buf, 1, NULL, PETSC_CHAR));
64 }
65 /* Read vertex information */
66 PetscCall(PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match));
67 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
68 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
69 PetscCall(PetscStrncmp(line, "vertex", PETSC_MAX_PATH_LEN, &match));
70 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
71 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
72 snum = sscanf(line, "%d", &Nv);
73 PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
74 match = PETSC_TRUE;
75 while (match) {
76 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
77 PetscCall(PetscStrncmp(line, "property", PETSC_MAX_PATH_LEN, &match));
78 if (match) {
79 PetscBool matchB;
80
81 PetscCall(PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING));
82 PetscCheck(Nvp < 16, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot handle more than 16 property statements in PLY file header: %s", line);
83 snum = sscanf(line, "%s %s", ntype, name);
84 PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
85 PetscCall(PetscStrncmp(ntype, "float32", 16, &matchB));
86 if (matchB) {
87 vtype[Nvp] = 'f';
88 } else {
89 PetscCall(PetscStrncmp(ntype, "int32", 16, &matchB));
90 if (matchB) {
91 vtype[Nvp] = 'd';
92 } else {
93 PetscCall(PetscStrncmp(ntype, "uint8", 16, &matchB));
94 PetscCheck(matchB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse type in PLY file header: %s", line);
95 vtype[Nvp] = 'c';
96 }
97 }
98 PetscCall(PetscStrncmp(name, "x", 16, &matchB));
99 if (matchB) xi = Nvp;
100 PetscCall(PetscStrncmp(name, "y", 16, &matchB));
101 if (matchB) yi = Nvp;
102 PetscCall(PetscStrncmp(name, "z", 16, &matchB));
103 if (matchB) zi = Nvp;
104 ++Nvp;
105 }
106 }
107 /* Read cell information */
108 PetscCall(PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match));
109 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
110 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
111 PetscCall(PetscStrncmp(line, "face", PETSC_MAX_PATH_LEN, &match));
112 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
113 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
114 snum = sscanf(line, "%d", &Nc);
115 PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
116 PetscCall(PetscViewerRead(viewer, line, 5, NULL, PETSC_STRING));
117 snum = sscanf(line, "property list %s %s %s", ntype, itype, name);
118 PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line);
119 PetscCall(PetscStrncmp(ntype, "uint8", 1024, &match));
120 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid size type in PLY file header: %s", line);
121 PetscCall(PetscStrncmp(name, "vertex_indices", 1024, &match));
122 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid property in PLY file header: %s", line);
123 /* Header should terminate */
124 PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
125 PetscCall(PetscStrncmp(line, "end_header", PETSC_MAX_PATH_LEN, &match));
126 PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file: %s", line);
127 } else {
128 Nc = Nv = 0;
129 }
130 PetscCall(DMCreate(comm, dm));
131 PetscCall(DMSetType(*dm, DMPLEX));
132 PetscCall(DMPlexSetChart(*dm, 0, Nc + Nv));
133 PetscCall(DMSetDimension(*dm, dim));
134 PetscCall(DMSetCoordinateDim(*dm, cdim));
135 /* Read coordinates */
136 PetscCall(DMGetCoordinateSection(*dm, &coordSection));
137 PetscCall(PetscSectionSetNumFields(coordSection, 1));
138 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, cdim));
139 PetscCall(PetscSectionSetChart(coordSection, Nc, Nc + Nv));
140 for (v = Nc; v < Nc + Nv; ++v) {
141 PetscCall(PetscSectionSetDof(coordSection, v, cdim));
142 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, cdim));
143 }
144 PetscCall(PetscSectionSetUp(coordSection));
145 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
146 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
147 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
148 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
149 PetscCall(VecSetBlockSize(coordinates, cdim));
150 PetscCall(VecSetType(coordinates, VECSTANDARD));
151 PetscCall(VecGetArray(coordinates, &coords));
152 if (rank == 0) {
153 float rbuf[1];
154 int ibuf[1];
155
156 for (v = 0; v < Nv; ++v) {
157 for (p = 0; p < Nvp; ++p) {
158 if (vtype[p] == 'f') {
159 PetscCall(PetscViewerRead(viewer, &rbuf, 1, NULL, PETSC_FLOAT));
160 if (byteSwap) PetscCall(PetscByteSwap(&rbuf, PETSC_FLOAT, 1));
161 if (p == xi) coords[v * cdim + 0] = rbuf[0];
162 else if (p == yi) coords[v * cdim + 1] = rbuf[0];
163 else if (p == zi) coords[v * cdim + 2] = rbuf[0];
164 } else if (vtype[p] == 'd') {
165 PetscCall(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_INT));
166 if (byteSwap) PetscCall(PetscByteSwap(&ibuf, PETSC_INT, 1));
167 } else if (vtype[p] == 'c') {
168 PetscCall(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR));
169 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid vertex property type in PLY file");
170 }
171 }
172 }
173 PetscCall(VecRestoreArray(coordinates, &coords));
174 PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
175 PetscCall(VecDestroy(&coordinates));
176 /* Read topology */
177 if (rank == 0) {
178 char ibuf[1];
179 PetscInt vbuf[16], corners;
180
181 /* Assume same cells */
182 PetscCall(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR));
183 corners = ibuf[0];
184 for (c = 0; c < Nc; ++c) PetscCall(DMPlexSetConeSize(*dm, c, corners));
185 PetscCall(DMSetUp(*dm));
186 for (c = 0; c < Nc; ++c) {
187 if (c > 0) PetscCall(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR));
188 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);
189 PetscCall(PetscViewerRead(viewer, &vbuf, ibuf[0], NULL, PETSC_INT));
190 if (byteSwap) PetscCall(PetscByteSwap(&vbuf, PETSC_INT, ibuf[0]));
191 for (v = 0; v < ibuf[0]; ++v) vbuf[v] += Nc;
192 PetscCall(DMPlexSetCone(*dm, c, vbuf));
193 }
194 }
195 PetscCall(DMPlexSymmetrize(*dm));
196 PetscCall(DMPlexStratify(*dm));
197 PetscCall(PetscViewerDestroy(&viewer));
198 if (interpolate) {
199 DM idm;
200
201 PetscCall(DMPlexInterpolate(*dm, &idm));
202 PetscCall(DMDestroy(dm));
203 *dm = idm;
204 }
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207