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 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 34 CHKERRQ(PetscViewerCreate(comm, &viewer)); 35 CHKERRQ(PetscViewerSetType(viewer, PETSCVIEWERBINARY)); 36 CHKERRQ(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 37 CHKERRQ(PetscViewerFileSetName(viewer, filename)); 38 if (rank == 0) { 39 PetscBool isAscii, isBinaryBig, isBinaryLittle; 40 41 /* Check for PLY file */ 42 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 43 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 47 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 50 CHKERRQ(PetscStrncmp(line, "ascii", PETSC_MAX_PATH_LEN, &isAscii)); 51 CHKERRQ(PetscStrncmp(line, "binary_big_endian", PETSC_MAX_PATH_LEN, &isBinaryBig)); 52 CHKERRQ(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 else if (isBinaryLittle) byteSwap = PETSC_TRUE; 55 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 56 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 63 CHKERRQ(PetscStrncmp(line, "comment", PETSC_MAX_PATH_LEN, &match)); 64 if (match) while (buf != '\n') CHKERRQ(PetscViewerRead(viewer, &buf, 1, NULL, PETSC_CHAR)); 65 } 66 /* Read vertex information */ 67 CHKERRQ(PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match)); 68 PetscCheck(match,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 69 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 70 CHKERRQ(PetscStrncmp(line, "vertex", PETSC_MAX_PATH_LEN, &match)); 71 PetscCheck(match,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 72 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 73 snum = sscanf(line, "%d", &Nv); 74 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 75 match = PETSC_TRUE; 76 while (match) { 77 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 78 CHKERRQ(PetscStrncmp(line, "property", PETSC_MAX_PATH_LEN, &match)); 79 if (match) { 80 PetscBool matchB; 81 82 CHKERRQ(PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING)); 83 PetscCheckFalse(Nvp >= 16,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot handle more than 16 property statements in PLY file header: %s", line); 84 snum = sscanf(line, "%s %s", ntype, name); 85 PetscCheckFalse(snum != 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 86 CHKERRQ(PetscStrncmp(ntype, "float32", 16, &matchB)); 87 if (matchB) { 88 vtype[Nvp] = 'f'; 89 } else { 90 CHKERRQ(PetscStrncmp(ntype, "int32", 16, &matchB)); 91 if (matchB) { 92 vtype[Nvp] = 'd'; 93 } else { 94 CHKERRQ(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 CHKERRQ(PetscStrncmp(name, "x", 16, &matchB)); 101 if (matchB) {xi = Nvp;} 102 CHKERRQ(PetscStrncmp(name, "y", 16, &matchB)); 103 if (matchB) {yi = Nvp;} 104 CHKERRQ(PetscStrncmp(name, "z", 16, &matchB)); 105 if (matchB) {zi = Nvp;} 106 ++Nvp; 107 } 108 } 109 /* Read cell information */ 110 CHKERRQ(PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match)); 111 PetscCheck(match,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 112 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 113 CHKERRQ(PetscStrncmp(line, "face", PETSC_MAX_PATH_LEN, &match)); 114 PetscCheck(match,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 115 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 116 snum = sscanf(line, "%d", &Nc); 117 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 118 CHKERRQ(PetscViewerRead(viewer, line, 5, NULL, PETSC_STRING)); 119 snum = sscanf(line, "property list %s %s %s", ntype, itype, name); 120 PetscCheckFalse(snum != 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 121 CHKERRQ(PetscStrncmp(ntype, "uint8", 1024, &match)); 122 PetscCheck(match,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid size type in PLY file header: %s", line); 123 CHKERRQ(PetscStrncmp(name, "vertex_indices", 1024, &match)); 124 PetscCheck(match,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid property in PLY file header: %s", line); 125 /* Header should terminate */ 126 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 127 CHKERRQ(PetscStrncmp(line, "end_header", PETSC_MAX_PATH_LEN, &match)); 128 PetscCheck(match,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file: %s", line); 129 } else { 130 Nc = Nv = 0; 131 } 132 CHKERRQ(DMCreate(comm, dm)); 133 CHKERRQ(DMSetType(*dm, DMPLEX)); 134 CHKERRQ(DMPlexSetChart(*dm, 0, Nc+Nv)); 135 CHKERRQ(DMSetDimension(*dm, dim)); 136 CHKERRQ(DMSetCoordinateDim(*dm, cdim)); 137 /* Read coordinates */ 138 CHKERRQ(DMGetCoordinateSection(*dm, &coordSection)); 139 CHKERRQ(PetscSectionSetNumFields(coordSection, 1)); 140 CHKERRQ(PetscSectionSetFieldComponents(coordSection, 0, cdim)); 141 CHKERRQ(PetscSectionSetChart(coordSection, Nc, Nc + Nv)); 142 for (v = Nc; v < Nc+Nv; ++v) { 143 CHKERRQ(PetscSectionSetDof(coordSection, v, cdim)); 144 CHKERRQ(PetscSectionSetFieldDof(coordSection, v, 0, cdim)); 145 } 146 CHKERRQ(PetscSectionSetUp(coordSection)); 147 CHKERRQ(PetscSectionGetStorageSize(coordSection, &coordSize)); 148 CHKERRQ(VecCreate(PETSC_COMM_SELF, &coordinates)); 149 CHKERRQ(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 150 CHKERRQ(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 151 CHKERRQ(VecSetBlockSize(coordinates, cdim)); 152 CHKERRQ(VecSetType(coordinates, VECSTANDARD)); 153 CHKERRQ(VecGetArray(coordinates, &coords)); 154 if (rank == 0) { 155 float rbuf[1]; 156 int ibuf[1]; 157 158 for (v = 0; v < Nv; ++v) { 159 for (p = 0; p < Nvp; ++p) { 160 if (vtype[p] == 'f') { 161 CHKERRQ(PetscViewerRead(viewer, &rbuf, 1, NULL, PETSC_FLOAT)); 162 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_INT)); 168 if (byteSwap) CHKERRQ(PetscByteSwap(&ibuf, PETSC_INT, 1)); 169 } else if (vtype[p] == 'c') { 170 CHKERRQ(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 CHKERRQ(VecRestoreArray(coordinates, &coords)); 176 CHKERRQ(DMSetCoordinatesLocal(*dm, coordinates)); 177 CHKERRQ(VecDestroy(&coordinates)); 178 /* Read topology */ 179 if (rank == 0) { 180 char ibuf[1]; 181 PetscInt vbuf[16], corners; 182 183 /* Assume same cells */ 184 CHKERRQ(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR)); 185 corners = ibuf[0]; 186 for (c = 0; c < Nc; ++c) CHKERRQ(DMPlexSetConeSize(*dm, c, corners)); 187 CHKERRQ(DMSetUp(*dm)); 188 for (c = 0; c < Nc; ++c) { 189 if (c > 0) { 190 CHKERRQ(PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR)); 191 } 192 PetscCheckFalse(ibuf[0] != corners,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "All cells must have the same number of vertices in PLY file: %D != %D", ibuf[0], corners); 193 CHKERRQ(PetscViewerRead(viewer, &vbuf, ibuf[0], NULL, PETSC_INT)); 194 if (byteSwap) CHKERRQ(PetscByteSwap(&vbuf, PETSC_INT, ibuf[0])); 195 for (v = 0; v < ibuf[0]; ++v) vbuf[v] += Nc; 196 CHKERRQ(DMPlexSetCone(*dm, c, vbuf)); 197 } 198 } 199 CHKERRQ(DMPlexSymmetrize(*dm)); 200 CHKERRQ(DMPlexStratify(*dm)); 201 CHKERRQ(PetscViewerDestroy(&viewer)); 202 if (interpolate) { 203 DM idm; 204 205 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 206 CHKERRQ(DMDestroy(dm)); 207 *dm = idm; 208 } 209 PetscFunctionReturn(0); 210 } 211