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: 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)
64: while (buf != '\n') PetscViewerRead(viewer, &buf, 1, NULL, PETSC_CHAR);
65: }
66: /* Read vertex information */
67: PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match);
69: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
70: PetscStrncmp(line, "vertex", PETSC_MAX_PATH_LEN, &match);
72: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
73: snum = sscanf(line, "%d", &Nv);
75: match = PETSC_TRUE;
76: while (match) {
77: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
78: PetscStrncmp(line, "property", PETSC_MAX_PATH_LEN, &match);
79: if (match) {
80: PetscBool matchB;
82: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
84: snum = sscanf(line, "%s %s", ntype, name);
86: PetscStrncmp(ntype, "float32", 16, &matchB);
87: if (matchB) {
88: vtype[Nvp] = 'f';
89: } else {
90: PetscStrncmp(ntype, "int32", 16, &matchB);
91: if (matchB) {
92: vtype[Nvp] = 'd';
93: } else {
94: PetscStrncmp(ntype, "uint8", 16, &matchB);
95: if (matchB) {
96: vtype[Nvp] = 'c';
97: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse type in PLY file header: %s", line);
98: }
99: }
100: PetscStrncmp(name, "x", 16, &matchB);
101: if (matchB) xi = Nvp;
102: PetscStrncmp(name, "y", 16, &matchB);
103: if (matchB) yi = Nvp;
104: PetscStrncmp(name, "z", 16, &matchB);
105: if (matchB) zi = Nvp;
106: ++Nvp;
107: }
108: }
109: /* Read cell information */
110: PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match);
112: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
113: PetscStrncmp(line, "face", PETSC_MAX_PATH_LEN, &match);
115: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
116: snum = sscanf(line, "%d", &Nc);
118: PetscViewerRead(viewer, line, 5, NULL, PETSC_STRING);
119: snum = sscanf(line, "property list %s %s %s", ntype, itype, name);
121: PetscStrncmp(ntype, "uint8", 1024, &match);
123: PetscStrncmp(name, "vertex_indices", 1024, &match);
125: /* Header should terminate */
126: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
127: PetscStrncmp(line, "end_header", PETSC_MAX_PATH_LEN, &match);
129: } else {
130: Nc = Nv = 0;
131: }
132: DMCreate(comm, dm);
133: DMSetType(*dm, DMPLEX);
134: DMPlexSetChart(*dm, 0, Nc + Nv);
135: DMSetDimension(*dm, dim);
136: DMSetCoordinateDim(*dm, cdim);
137: /* Read coordinates */
138: DMGetCoordinateSection(*dm, &coordSection);
139: PetscSectionSetNumFields(coordSection, 1);
140: PetscSectionSetFieldComponents(coordSection, 0, cdim);
141: PetscSectionSetChart(coordSection, Nc, Nc + Nv);
142: for (v = Nc; v < Nc + Nv; ++v) {
143: PetscSectionSetDof(coordSection, v, cdim);
144: PetscSectionSetFieldDof(coordSection, v, 0, cdim);
145: }
146: PetscSectionSetUp(coordSection);
147: PetscSectionGetStorageSize(coordSection, &coordSize);
148: VecCreate(PETSC_COMM_SELF, &coordinates);
149: PetscObjectSetName((PetscObject)coordinates, "coordinates");
150: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
151: VecSetBlockSize(coordinates, cdim);
152: VecSetType(coordinates, VECSTANDARD);
153: VecGetArray(coordinates, &coords);
154: if (rank == 0) {
155: float rbuf[1];
156: int ibuf[1];
158: for (v = 0; v < Nv; ++v) {
159: for (p = 0; p < Nvp; ++p) {
160: if (vtype[p] == 'f') {
161: PetscViewerRead(viewer, &rbuf, 1, NULL, PETSC_FLOAT);
162: if (byteSwap) PetscByteSwap(&rbuf, PETSC_FLOAT, 1);
163: if (p == xi) coords[v * cdim + 0] = rbuf[0];
164: else if (p == yi) coords[v * cdim + 1] = rbuf[0];
165: else if (p == zi) coords[v * cdim + 2] = rbuf[0];
166: } else if (vtype[p] == 'd') {
167: PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_INT);
168: if (byteSwap) PetscByteSwap(&ibuf, PETSC_INT, 1);
169: } else if (vtype[p] == 'c') {
170: PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);
171: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid vertex property type in PLY file");
172: }
173: }
174: }
175: VecRestoreArray(coordinates, &coords);
176: DMSetCoordinatesLocal(*dm, coordinates);
177: VecDestroy(&coordinates);
178: /* Read topology */
179: if (rank == 0) {
180: char ibuf[1];
181: PetscInt vbuf[16], corners;
183: /* Assume same cells */
184: PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);
185: corners = ibuf[0];
186: for (c = 0; c < Nc; ++c) DMPlexSetConeSize(*dm, c, corners);
187: DMSetUp(*dm);
188: for (c = 0; c < Nc; ++c) {
189: if (c > 0) PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);
191: PetscViewerRead(viewer, &vbuf, ibuf[0], NULL, PETSC_INT);
192: if (byteSwap) PetscByteSwap(&vbuf, PETSC_INT, ibuf[0]);
193: for (v = 0; v < ibuf[0]; ++v) vbuf[v] += Nc;
194: DMPlexSetCone(*dm, c, vbuf);
195: }
196: }
197: DMPlexSymmetrize(*dm);
198: DMPlexStratify(*dm);
199: PetscViewerDestroy(&viewer);
200: if (interpolate) {
201: DM idm;
203: DMPlexInterpolate(*dm, &idm);
204: DMDestroy(dm);
205: *dm = idm;
206: }
207: return 0;
208: }