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