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 @*/ 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