1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 /*@C 5 DMPlexCreatePLYFromFile - Create a `DMPLEX` mesh from a PLY file. 6 7 Input Parameters: 8 + comm - The MPI communicator 9 . filename - Name of the .med file 10 - interpolate - Create faces and edges in the mesh 11 12 Output Parameter: 13 . dm - The DM object representing the mesh 14 15 References: 16 . * - https://en.wikipedia.org/wiki/PLY_(file_format) 17 18 Level: beginner 19 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; 33 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; 42 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; 84 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]; 160 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; 185 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; 205 206 PetscCall(DMPlexInterpolate(*dm, &idm)); 207 PetscCall(DMDestroy(dm)); 208 *dm = idm; 209 } 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212